ebcs1.F File Reference

#include "implicit_f.inc"
#include "com08_c.inc"
#include "scr11_c.inc"

Go to the source code of this file.

Functions/Subroutines
subroutine	ebcs1 (nseg, iseg, segvar, a, v, x, liste, nod, irect, vo, po, la, fv, ms, stifn, ebcs_tab, iebcs)

Function/Subroutine Documentation

ebcs1()

subroutine ebcs1	(	integer	nseg,
		integer, dimension(nseg)	iseg,
		type(t_segvar)	segvar,
			a,
			v,
			x,
		integer, dimension(nod)	liste,
		integer	nod,
		integer, dimension(4,nseg)	irect,
			vo,
			po,
			la,
			fv,
			ms,
			stifn,
		type(t_ebcs_tab), intent(in), target	ebcs_tab,
		integer, intent(in)	iebcs )

Definition at line 31 of file ebcs1.F.

      USE ebcs_mod
      USE segvar_mod
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "com08_c.inc"
#include      "scr11_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER NSEG,NOD,ISEG(NSEG),LISTE(NOD),IRECT(4,NSEG)
      my_real
     .        a(3,*),x(3,*),v(3,*),la(3,nod),
     .        vo(nod),po(nod),fv(*),ms(*),stifn(*)
      TYPE(t_ebcs_tab), TARGET, INTENT(IN) :: EBCS_TAB
      INTEGER, INTENT(IN) :: IEBCS
      TYPE(t_segvar) :: SEGVAR
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER IPRES,IRHO,I,IS,KSEG,N1,N2,N3,N4,NG1,NG2,NG3,NG4,
     .        N,IENER
      my_real
     .       orient,pres,rho,c,lcar,r1,r2,roc,alp,fac,
     .       x13,y13,z13,x24,y24,z24,nx,ny,nz,s,
     .       roou,enou,vmx,vmy,vmz,fluxi,fluxo,vn,pn,du,dp,p,
     .       ener,dpdv
      CLASS(t_ebcs), POINTER :: EBCS
C-----------------------------------------------
C Presssion et densite imposes
      ebcs => ebcs_tab%tab(iebcs)%poly
      c = zero
      lcar = zero
      ener = zero
      pres = zero
      rho = zero
      r1 = zero
      r2 = zero
      ipres = 0 
      irho = 0 
      iener = 0 
      SELECT TYPE (ebcs)
      TYPE IS (t_ebcs_pres)
        ipres=ebcs%ipres
        irho=ebcs%irho
        iener=ebcs%iener
        IF(ipres>0)THEN
           pres=ebcs%pres*fv(ipres)
        ELSE
           pres=ebcs%pres
        ENDIF
        IF(irho>0)THEN
           rho=ebcs%rho*fv(irho)
        ELSE
           rho=ebcs%rho
        ENDIF
        IF(iener>0)THEN
           ener=ebcs%ener*fv(iener)
        ELSE
           ener=ebcs%ener
        ENDIF
        c=ebcs%c
        lcar=ebcs%lcar
        r1=ebcs%r1
        r2=ebcs%r2
      TYPE IS (t_ebcs_valvin)
        ipres=ebcs%ipres
        irho=ebcs%irho
        iener=ebcs%iener
        IF(ipres>0)THEN
           pres=ebcs%pres*fv(ipres)
        ELSE
           pres=ebcs%pres
        ENDIF
        IF(irho>0)THEN
           rho=ebcs%rho*fv(irho)
        ELSE
           rho=ebcs%rho
        ENDIF
        IF(iener>0)THEN
           ener=ebcs%ener*fv(iener)
        ELSE
           ener=ebcs%ener
        ENDIF
        c=ebcs%c
        lcar=ebcs%lcar
        r1=ebcs%r1
        r2=ebcs%r2
      TYPE IS (t_ebcs_valvout)
        ipres=ebcs%ipres
        irho=ebcs%irho
        iener=ebcs%iener
        IF(ipres>0)THEN
           pres=ebcs%pres*fv(ipres)
        ELSE
           pres=ebcs%pres
        ENDIF
        IF(irho>0)THEN
           rho=ebcs%rho*fv(irho)
        ELSE
           rho=ebcs%rho
        ENDIF
        IF(iener>0)THEN
           ener=ebcs%ener*fv(iener)
        ELSE
           ener=ebcs%ener
        ENDIF
        c=ebcs%c
        lcar=ebcs%lcar
        r1=ebcs%r1
        r2=ebcs%r2
      END SELECT
 
      roc=rho*c
c
c      write(6,*)R1,R2,LCAR
c
      alp=zero
      IF (lcar>0)alp=c*dt1/lcar
C SURFACE NORMALE NODALES
      DO i=1,nod
         la(1,i)=zero
         la(2,i)=zero
         la(3,i)=zero
      ENDDO
      DO is=1,nseg
        kseg=abs(iseg(is))
        orient=float(iseg(is)/kseg)
        n1=irect(1,is)
        n2=irect(2,is)
        n3=irect(3,is)
        n4=irect(4,is)
        IF(n4==0 .OR. n4==n3) THEN
          fac=one_over_6*orient
          n4=n3
        ELSE
          fac=one_over_8*orient
        ENDIF
c
        ng1=liste(n1)
        ng2=liste(n2)
        ng3=liste(n3)
        ng4=liste(n4)
        x13=x(1,ng3)-x(1,ng1)
        y13=x(2,ng3)-x(2,ng1)
        z13=x(3,ng3)-x(3,ng1)
        x24=x(1,ng4)-x(1,ng2)
        y24=x(2,ng4)-x(2,ng2)
        z24=x(3,ng4)-x(3,ng2)
c
        nx=(y13*z24-z13*y24)*fac
        ny=(z13*x24-x13*z24)*fac
        nz=(x13*y24-y13*x24)*fac
c
        la(1,n1)=la(1,n1)+nx
        la(2,n1)=la(2,n1)+ny
        la(3,n1)=la(3,n1)+nz
        la(1,n2)=la(1,n2)+nx
        la(2,n2)=la(2,n2)+ny
        la(3,n2)=la(3,n2)+nz
        la(1,n3)=la(1,n3)+nx
        la(2,n3)=la(2,n3)+ny
        la(3,n3)=la(3,n3)+nz 
C
        vmx=v(1,ng1)+v(1,ng2)+v(1,ng3)
        vmy=v(2,ng1)+v(2,ng2)+v(2,ng3)
        vmz=v(3,ng1)+v(3,ng2)+v(3,ng3)
        IF(n4/=n3) THEN
           la(1,n4)=la(1,n4)+nx
           la(2,n4)=la(2,n4)+ny
           la(3,n4)=la(3,n4)+nz
           vmx=vmx+v(1,ng4)
           vmy=vmy+v(2,ng4)
           vmz=vmz+v(3,ng4)
        ENDIF
C
c bilan masse et energie totale
c
        roou = segvar%RHO(kseg)
        enou = segvar%EINT(kseg)
C
        fluxo=(vmx*nx+vmy*ny+vmz*nz)*dt1
        fluxi=min(fluxo,zero)
        fluxo=max(fluxo,zero)
        dmf=dmf-fluxo*roou-fluxi*rho
        def=def-fluxo*enou-fluxi*ener
C
C stockage densite et energie entrante dans buffer facette
C
        segvar%RHO(kseg)=rho
        segvar%EINT(kseg)=ener
      ENDDO
C
      DO i=1,nod
        n=liste(i)
        s=sqrt(la(1,i)**2+la(2,i)**2+la(3,i)**2)
        vn=(v(1,n)*la(1,i)+v(2,n)*la(2,i)+v(3,n)*la(3,i))/s
C ajout resistance
        pn=pres+r1*vn+r2*vn*abs(vn)
        dpdv=roc+r1+two*r2*abs(vn)
c condition darret
c        write(6,*)i,n,pn,pres,-0.5*RHO*VN**2
        IF(vn<0)THEN
          pn=pn-half*rho*vn**2
          dpdv=dpdv-rho*vn
        ENDIF
c frontiere silencieuse
        IF(tt>0)THEN
          du=roc*(vn-vo(i))
        ELSE
          du=zero
          po(i)=pn
        ENDIF
        dp=alp*(pn-po(i))
c        write(6,*)'vitesse',vn,vo(i)
        vo(i)=vn
        p=po(i)+dp+du
        IF(c==zero)p=pn
c        write(6,*)'ebcs1 pression',p,'vn',vn,'dp',dp,'du',du
c
        a(1,n)=a(1,n)-p*la(1,i)
        a(2,n)=a(2,n)-p*la(2,i)
        a(3,n)=a(3,n)-p*la(3,i)
C contribution des forces a la perte d'energie globale
        def=def-half*(po(i)+p)*dt1*vn*s
        po(i)=p
        stifn(n)=stifn(n)+(two*(s*dpdv)**2)/ms(n)
c        write(6,*)'stifn',STIFN(N)
      ENDDO
c
      RETURN