46
47
48
49 USE ebcs_mod
51 USE output_mod , ONLY : output_
52
53
54
55#include "implicit_f.inc"
56
57
58
59
60
61
62 INTEGER NSEG,NOD,ISEG(NSEG),LISTE(NOD),IRECT(4,NSEG)
63 my_real :: a(3,*),x(3,*),v(3,*),la(3,nod),vo(nod),po(nod),fv(*),ms(*),stifn(*)
64 TYPE(t_ebcs_tab), TARGET, INTENT(IN) :: EBCS_TAB
65 INTEGER, INTENT(IN) :: IEBCS
66 TYPE(t_segvar) :: SEGVAR
67 TYPE(OUTPUT_), INTENT(INOUT) :: OUTPUT
70
71
72
73 INTEGER :: IPRES,IRHO,I,IS,KSEG,N1,N2,N3,N4,NG1,NG2,NG3,NG4,N,IENER
74 my_real :: orient,pres,rho,c,lcar,r1,r2,roc,alp,fac,
75 . x13,y13,z13,x24,y24,z24,nx,ny,nz,s,
76 . roou,enou,vmx,vmy,vmz,fluxi,fluxo,vn,pn,du,dp,p,
77 . ener,dpdv
78 my_real :: de, de_in, de_out, dm_in, dm_out
79 CLASS(t_ebcs), POINTER :: EBCS
80
81
82 ebcs => ebcs_tab%tab(iebcs)%poly
83 c = zero
84 lcar = zero
85 ener = zero
86 pres = zero
87 rho = zero
88 r1 = zero
89 r2 = zero
90 ipres = 0
91 irho = 0
92 iener = 0
93 de_in = zero
94 de_out = zero
95 dm_in = zero
96 dm_out = zero
97
98 SELECT TYPE (ebcs)
99 TYPE IS (t_ebcs_pres)
100 ipres=ebcs%ipres
101 irho=ebcs%irho
102 iener=ebcs%iener
103 IF(ipres > 0)THEN
104 pres=ebcs%pres*fv(ipres)
105 ELSE
106 pres=ebcs%pres
107 ENDIF
108 IF(irho > 0)THEN
109 rho=ebcs%rho*fv(irho)
110 ELSE
111 rho=ebcs%rho
112 ENDIF
113 IF(iener > 0)THEN
114 ener=ebcs%ener*fv(iener)
115 ELSE
116 ener=ebcs%ener
117 ENDIF
118 c=ebcs%c
119 lcar=ebcs%lcar
120 r1=ebcs%r1
121 r2=ebcs%r2
122 TYPE IS (t_ebcs_valvin)
123 ipres=ebcs%ipres
124 irho=ebcs%irho
125 iener=ebcs%iener
126 IF(ipres > 0)THEN
127 pres=ebcs%pres*fv(ipres)
128 ELSE
129 pres=ebcs%pres
130 ENDIF
131 IF(irho > 0)THEN
132 rho=ebcs%rho*fv(irho)
133 ELSE
134 rho=ebcs%rho
135 ENDIF
136 IF(iener > 0)THEN
137 ener=ebcs%ener*fv(iener)
138 ELSE
139 ener=ebcs%ener
140 ENDIF
141 c=ebcs%c
142 lcar=ebcs%lcar
143 r1=ebcs%r1
144 r2=ebcs%r2
145 TYPE IS (t_ebcs_valvout)
146 ipres=ebcs%ipres
147 irho=ebcs%irho
148 iener=ebcs%iener
149 IF(ipres > 0)THEN
150 pres=ebcs%pres*fv(ipres)
151 ELSE
152 pres=ebcs%pres
153 ENDIF
154 IF(irho > 0)THEN
155 rho=ebcs%rho*fv(irho)
156 ELSE
157 rho=ebcs%rho
158 ENDIF
159 IF(iener > 0)THEN
160 ener=ebcs%ener*fv(iener)
161 ELSE
162 ener=ebcs%ener
163 ENDIF
164 c=ebcs%c
165 lcar=ebcs%lcar
166 r1=ebcs%r1
167 r2=ebcs%r2
168 END SELECT
169
170 roc=rho*c
171
172 alp=zero
173 IF (lcar > zero)alp=c*dt1/lcar
174
175 DO i=1,nod
176 la(1,i)=zero
177 la(2,i)=zero
178 la(3,i)=zero
179 ENDDO
180
181 DO is=1,nseg
182 kseg=abs(iseg(is))
183 orient=float(iseg(is)/kseg)
184 n1=irect(1,is)
185 n2=irect(2,is)
186 n3=irect(3,is)
187 n4=irect(4,is)
188 IF(n4==0 .OR. n4==n3) THEN
189 fac=one_over_6*orient
190 n4=n3
191 ELSE
192 fac=one_over_8*orient
193 ENDIF
194
195 ng1=liste(n1)
196 ng2=liste(n2)
197 ng3=liste(n3)
198 ng4=liste(n4)
199 x13=x(1,ng3)-x(1,ng1)
200 y13=x(2,ng3)-x(2,ng1)
201 z13=x(3,ng3)-x(3,ng1)
202 x24=x(1,ng4)-x(1,ng2)
203 y24=x(2,ng4)-x(2,ng2)
204 z24=x(3,ng4)-x(3,ng2)
205
206 nx=(y13*z24-z13*y24)*fac
207 ny=(z13*x24-x13*z24)*fac
208 nz=(x13*y24-y13*x24)*fac
209
210 la(1,n1)=la(1,n1)+nx
211 la(2,n1)=la(2,n1)+ny
212 la(3,n1)=la(3,n1)+nz
213 la(1,n2)=la(1,n2)+nx
214 la(2,n2)=la(2,n2)+ny
215 la(3,n2)=la(3,n2)+nz
216 la(1,n3)=la(1,n3)+nx
217 la(2,n3)=la(2,n3)+ny
218 la(3,n3)=la(3,n3)+nz
219
220 vmx=v(1,ng1)+v(1,ng2)+v(1,ng3)
221 vmy=v(2,ng1)+v(2,ng2)+v(2,ng3)
222 vmz=v(3,ng1)+v(3,ng2)+v(3,ng3)
223 IF(n4 /= n3) THEN
224 la(1,n4)=la(1,n4)+nx
225 la(2,n4)=la(2,n4)+ny
226 la(3,n4)=la(3,n4)+nz
227 vmx=vmx+v(1,ng4)
228 vmy=vmy+v(2,ng4)
229 vmz=vmz+v(3,ng4)
230 ENDIF
231
232
233
234 roou = segvar%RHO(kseg)
235 enou = segvar%EINT(kseg)
236
237 fluxo=(vmx*nx+vmy*ny+vmz*nz)*dt1
238 fluxi=
min(fluxo,zero)
239 fluxo=
max(fluxo,zero)
240 dm_out=dm_out-fluxo*roou
241 dm_in=dm_in-fluxi*rho
242 de_out=de_out-fluxo*enou
243 de_in=de_in-fluxi*ener
244
245
246
247 segvar%RHO(kseg)=rho
248 segvar%EINT(kseg)=ener
249 ENDDO
250
251
252 DO i=1,nod
253 n=liste(i)
254 s=sqrt(la(1,i)**2+la(2,i)**2+la(3,i)**2)
255 vn=(v(1,n)*la(1,i)+v(2,n)*la(2,i)+v(3,n)*la(3,i))/s
256
257 pn=pres+r1*vn+r2*vn*abs(vn)
258 dpdv=roc+r1+two*r2*abs(vn)
259
260
261 IF(vn<0)THEN
262 pn=pn-half*rho*vn**2
263 dpdv=dpdv-rho*vn
264 ENDIF
265
266 IF(time > zero)THEN
267 du=roc*(vn-vo(i))
268 ELSE
269 du=zero
270 po(i)=pn
271 ENDIF
272 dp=alp*(pn-po(i))
273
274 vo(i)=vn
275 p=po(i)+dp+du
276 IF(c == zero)p=pn
277
278
279 a(1,n)=a(1,n)-p*la(1,i)
280 a(2,n)=a(2,n)-p*la(2,i)
281 a(3,n)=a(3,n)-p*la(3,i)
282
283 de = -half*(po(i)+p)*dt1*vn*s
284 de_in = de_in +
max(de, zero)
285 de_out = de_out +
min(de, zero)
286 po(i)=p
287 stifn(n)=stifn(n)+(two*(s*dpdv)**2)/ms(n)
288
289 ENDDO
290
291
292 output%DATA%INOUT%DM_IN = output%DATA%INOUT%DM_IN + dm_in
293 output%DATA%INOUT%DM_OUT = output%DATA%INOUT%DM_OUT + dm_out
294 output%DATA%INOUT%DE_IN = output%DATA%INOUT%DE_IN + de_in
295 output%DATA%INOUT%DE_OUT = output%DATA%INOUT%DE_OUT + de_out
296
297
298 RETURN