32
33
34
35
36
37
38
39#include "implicit_f.inc"
40
41
42
43#include "mvsiz_p.inc"
44
45
46
47#include "vect01_c.inc"
48#include "param_c.inc"
49
50
51
52 my_real,
INTENT(IN) :: pm(npropm, *), geo(npropg, *), aire(*), vol(*)
53 my_real,
INTENT(INOUT) :: dtx(*)
54 my_real,
INTENT(IN),
DIMENSION(:),
TARGET :: bufmat(*), deltax(*)
55 INTEGER, INTENT(IN) :: PID(*),MAT(*),IPM(NPROPMI, *)
56
57
58
59 INTEGER :: I, MX, IADBUF,IFLG(MVSIZ)
60 my_real,
DIMENSION(:),
POINTER :: uparam
61
63 . ssp(mvsiz) , dpdm(mvsiz) , rho0(mvsiz) ,
64 . bulk , c1 , p ,
65 . mas(mvsiz), mu1p1(mvsiz),mu1p2(mvsiz),rho1(mvsiz),rho2(mvsiz),
66 . gam,rho10,rho20,p0,vfrac1,vfrac2,uvar1,a1,r1,visa1,visb1,visa2,visb2,
67 . b1,b2,mas1,mas2,ssp1,ssp2,vis(mvsiz),a
68
69 INTEGER ISOLVER
70
71
72
73
74
75
76
77 iadbuf = ipm(7,mat(1))
78 uparam =>bufmat(iadbuf:iadbuf+17)
79 isolver = uparam(17)
80 a1 = uparam(10)
81 rho10 = uparam(11)
82 rho20 = uparam(12)
83 r1 = uparam(6)
84 gam = uparam(5)
85 p0 = uparam(9)
86 c1 = uparam(4)
87 visa1 = uparam(1)
88 visb1 = uparam(3)
89 visa2 = uparam(13)
90 visb2 = uparam(15)
91
92 DO i=lft,llt
93 rho0(i)= rho10 * a1 + (one-a1)*rho20
94 mas(i)= rho0(i)*vol(i)
95 ENDDO
96
97 DO i=lft,llt
98 IF(gam*c1>=em30)THEN !if liquid and gas correctly defined
99 rho1(i) = rho10
100 rho2(i) = rho20
101 ENDIF
102 ENDDO
103
104 IF (isolver==2) THEN
105 DO i=lft,llt
106 mas1 = a1 * vol(i)
107 IF (a1 < em10)THEN
108 mas1=zero
109 vfrac1=zero
110 vfrac2=one
111 ENDIF
112 IF (one-a1 < em10)THEN
113 mas2=zero
114 rho1=mas(i)/vol(i)
115 vfrac1=one
116 vfrac2=zero
117 ENDIF
118 a = (rho0(i)-rho2(i))/(rho1(i)-rho2(i))
119 vfrac1=a
120 IF(vfrac1<em20)vfrac1=zero
121 vfrac2 = one-vfrac1
122 ssp1=c1
123 ssp2=gam*p0*(rho2(i)/rho20)**gam
124 uvar1=a*rho1(i)
125 IF(ssp1>zero)THEN
126 ssp1 = vfrac1 / ssp1
127 ELSE
128 ssp1=zero
129 ENDIF
130 IF(ssp2>zero)THEN
131 ssp2 = vfrac2 / ssp2
132 ELSE
133 ssp2=zero
134 ENDIF
135 ssp(i) = ssp1 + ssp2
136 ssp(i) = sqrt(one / ssp(i) / rho0(i))
137 b1=uvar1
138 b2=rho0(i)-b1
139 vis(i)=(b1*rho1(i)*visa1 + b2*rho2(i)*visa2)/rho0(i)
140
141 ENDDO
142 ELSE
143 DO i=lft,llt
144 ssp(i)=sqrt(r1)
145 vis(i)=zero
146 ENDDO
147 ENDIF
148
149
150
151
152
153 IF(jsph==0)THEN
154 CALL dtel(ssp,pm,geo,pid,mat, rho0, vis, deltax, aire, vol, dtx)
155 ELSE
156 CALL dtsph(ssp,pm,geo,pid,mat, rho0, vis, deltax, vol, dtx)
157 ENDIF
158
159 RETURN
subroutine dtel(ssp, pm, geo, pid, mat, rho0, vis, deltax, aire, vol, dtx)
subroutine dtsph(ssp, pm, geo, pid, mat, rho0, vis, deltax, vol, dtx)