41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
65 USE matparam_def_mod
68
69
70
71#include "implicit_f.inc"
72
73
74
75#include "units_c.inc"
76#include "com01_c.inc"
77#include "param_c.inc"
78#include "inter22.inc"
79
80
81
82 TYPE (UNIT_TYPE_)INTENT(IN) ::UNITAB
83 my_real,
DIMENSION(NPROPM) ,
INTENT(INOUT) :: pm
84 my_real,
DIMENSION(100) ,
INTENT(INOUT) :: parmat
85 my_real,
DIMENSION(MAXUPARAM) ,
INTENT(INOUT) :: uparam
86 INTEGER, DIMENSION(MAXFUNC) ,INTENT(INOUT) :: IFUNC
87 INTEGER, INTENT(INOUT) :: ISRATE,IMATVIS,NFUNC,MAXFUNC,MAXUPARAM,NUPARAM,NUVAR
88 TYPE(MLAW_TAG_),INTENT(INOUT) :: MTAG
89 INTEGER,INTENT(IN) :: ID
90 CHARACTER(LEN=NCHARTITLE) ,INTENT(IN) :: TITR
91 my_real,
INTENT(INOUT) :: stifint
92 TYPE(SUBMODEL_DATA),INTENT(IN) :: LSUBMODEL(NSUBMOD)
93 INTEGER,INTENT(IN) :: MAT_ID
94 TYPE(MATPARAM_STRUCT_) ,INTENT(INOUT) :: MATPARAM
95
96
97
98 LOGICAL :: IS_AVAILABLE,IS_ENCRYPTED
99 my_real c1,c2,p0,vis,visv,rho10,rho20,a1,pmin,gam,vis2,pshift,visv2,fac_m,fac_l,fac_t,fac_c,rho0,rhor
100 INTEGER I,ISOLVER,IRHO
101 character*63 chain1
102
103
104
105
106 israte = 0
107 is_encrypted = .false.
108 is_available = .false.
109
110
111
112
114 IF(trimat==0)trimat = -2
115
117
118 CALL hm_get_floatv(
'MAT_RHO' ,rho0 ,is_available, lsubmodel, unitab)
119 CALL hm_get_floatv(
'Refer_Rho' ,rhor ,is_available, lsubmodel, unitab)
120 CALL hm_get_floatv(
'MAT_PSH' ,pshift ,is_available, lsubmodel, unitab)
121
122
123 CALL hm_get_floatv(
'Lqud_Rho_l' ,rho10 ,is_available, lsubmodel, unitab)
124 CALL hm_get_floatv(
'C_l' ,c1 ,is_available, lsubmodel, unitab)
125 CALL hm_get_floatv(
'ALPHA1' ,a1 ,is_available, lsubmodel, unitab)
126 CALL hm_get_floatv(
'Nu_l' ,vis ,is_available, lsubmodel, unitab)
127 CALL hm_get_floatv(
'Bulk_Ratio_l' ,visv ,is_available, lsubmodel, unitab)
128
129
130 CALL hm_get_floatv('lqud_rho_g
' ,RHO20 ,IS_AVAILABLE, LSUBMODEL, UNITAB)
131 CALL HM_GET_FLOATV('lqud_gamma_bulk',GAM ,IS_AVAILABLE, LSUBMODEL, UNITAB)
132 CALL HM_GET_FLOATV('lqud_p0' ,P0 ,IS_AVAILABLE, LSUBMODEL, UNITAB)
133 CALL HM_GET_FLOATV('nu_g' ,VIS2 ,IS_AVAILABLE, LSUBMODEL, UNITAB)
134 CALL HM_GET_FLOATV('bulk_ratio_g' ,VISV2 ,IS_AVAILABLE, LSUBMODEL, UNITAB)
135
136 IRHO=0
137 IF(GAM*C1==0)IRHO=INT(P0)
138
139 !--------------------------------!
140 ! INPUT CHECK !
141 !--------------------------------!
142.OR. IF(A1<ZERO A1>ONE)THEN
143 chain1 = 'initial massic fraction must be between 0 and 1. '
144 CALL ANCMSG(MSGID=65,
145 . MSGTYPE=MSGERROR,
146 . ANMODE=ANINFO,
147 . I1=37,
148 . I2=ID,
149 . C1='error',
150 . C2=TITR,
151 . C3=chain1)
152 A1 = ONE
153 ENDIF
154
155.OR. IF(RHO10<ZERO RHO20<ZERO)THEN
156 chain1 = 'initial density cannot be negative. '
157 CALL ANCMSG(MSGID=65,
158 . MSGTYPE=MSGERROR,
159 . ANMODE=ANINFO,
160 . I1=37,
161 . I2=ID,
162 . C1='error',
163 . C2=TITR,
164 . C3=chain1)
165 IF(RHO10<ZERO)RHO10 = EM20
166 IF(RHO20<ZERO)RHO20 = EM20
167 ENDIF
168
169 !--------------------------------!
170 ! STORAGE !
171 !--------------------------------!
172
173 !===SOLVER VERSION
174 !ISOLVER = 1 : legacy solver
175 !ISOLVER = 2 : alternative solver (2018.0)
176 ISOLVER = 1
177 IF(INT22 > 0) ISOLVER = 2
178
179 !===NEW PRESSURE SHIFT PARAM (2018.0)
180 IF(PSHIFT==ZERO)THEN
181 PSHIFT=-P0
182 ELSE
183 PSHIFT = -PSHIFT
184 ENDIF
185
186 ! PSH = -P0 : relative pressure P_liq = C1.mu, P_gas = P0(1+mu)**gamma - P0
187 ! PSH = -1e-20 : total pressure P_liq = C0+C1.mu, P_gas = P0(1+mu)**gamma
188
189 !===UPARAM STORAGE
190 PMIN = - P0
191 RHO0 = RHO10 * A1 + (ONE-A1)*RHO20
192 RHOR = RHO0
193 PM(01) = RHO0
194 PM(89) = RHO0
195 PM(91) = RHO10
196 PM(31) = P0+PSHIFT
197 PM(104) = P0+PSHIFT
198
199 UPARAM(1) = VIS
200 UPARAM(2) = TWO*VIS
201 UPARAM(3) = (VISV-UPARAM(2))*THIRD
202 UPARAM(4) = C1
203 C2 = - PMIN * GAM
204 UPARAM(5) = GAM
205 UPARAM(6) = C1/RHO10
206 UPARAM(7) = GAM * P0 / RHO20
207 UPARAM(8) = PMIN
208 UPARAM(9) = P0
209 UPARAM(10) = A1
210 UPARAM(11) = RHO10
211 UPARAM(12) = RHO20
212 UPARAM(13) = VIS2
213 UPARAM(14) = TWO*VIS2
214 UPARAM(15) = (VISV2-UPARAM(14))*THIRD
215 UPARAM(16) = PSHIFT !pressure formulation : total or relative to P0
216 UPARAM(17) = ISOLVER
217
218 NUPARAM = 17
219 IFUNC(1) = IRHO
220 NFUNC = 1
221 NUVAR = 5
222 STIFINT = C1
223
224 ! MATPARAM keywords
225 CALL INIT_MAT_KEYWORD(MATPARAM,"HOOK")
226
227 ! Properties compatibility
228 CALL INIT_MAT_KEYWORD(MATPARAM,"SOLID_POROUS")
229
230 !--------------------------------!
231 ! PRINTOUT !
232 !--------------------------------!
233 IF(GAM*C1>0.)THEN
234 WRITE(IOUT,1011) TRIM(TITR),MAT_ID,37
235 WRITE(IOUT,1000)
236
237 IF(ISOLVER==2)WRITE(IOUT,1102)
238 IF(IS_ENCRYPTED)THEN
239 WRITE(IOUT,'(5x,a,//)')'confidential data'
240 ELSE
241 WRITE(IOUT,1012) RHO0
242 WRITE(IOUT,1100)RHO10,C1,A1,VIS,VISV,RHO20,GAM,P0,VIS2,VISV2,P0,-PSHIFT
243 ENDIF
244 ELSE
245 WRITE(IOUT,1001)
246 IF(IS_ENCRYPTED)THEN
247 WRITE(IOUT,'(5x,a,//)')'confidential data'
248 ELSE
249 WRITE(IOUT,1101)RHO10,RHO20,IRHO
250 END IF
251 END IF
252
253
254 RETURN
255 1000 FORMAT(
256 & 5X,' bi-phase liquid-gas law',/,
257 & 5X,' -----------------------',//)
258 1011 FORMAT(/
259 & 5X,A,/,
260 & 5X,'material number . . . . . . . . . . . .=',I10/,
261 & 5X,'material law. . . . . . . . . . . . . .=',I10/)
262 1012 FORMAT(
263 & 5X,'initial density . . . . . . . . . . . .=',1PG20.13/)
264 1001 FORMAT(
265 & 5X,' boundary
for bi-phase liquid-gas law
',/,
266 & 5X,' -------------------------------------',//)
267 1100 FORMAT(
268 & 5X,'liquid reference density. . . . . . . .=',1PG20.13/
269 & 5X,'liquid bulk stiffness . . . . . . . . .=',1PG20.13/
270 & 5X,'liquid mass ratio . . . . . . . . . . .=',1PG20.13/
271 & 5X,' (0. = full of gaz; 1. = full of liquid)'/
272 & 5X,'liquid shear kinematic viscosity . . . =',1PG20.13/
273 & 5X,'liquid bulk kinematic viscosity . . . =',1PG20.13//
274 & 5X,'gas reference density. . . . . . . .=',1PG20.13/
275 & 5X,'gas constant. . . . . . . . . . . . . .=',1PG20.13/
276 & 5X,'gas reference pressure. . . . . . . . .=',1PG20.13/
277 & 5X,'gas shear kinematic viscosity . . . . .=',1PG20.13/
278 & 5X,'gas bulk kinematic viscosity . . . . .=',1PG20.13//
279 & 5X,'reference pressure. . . . . . . . . . .=',1PG20.13/
280 & 5X,'pressure shift. . . . . . . . . . . . .=',1PG20.13//)
281 1101 FORMAT(
282 & 5X,'liquid reference density. . . . . . . .=',1PG20.13/
283 & 5X,'gas reference density. . . . . . . .=',1PG20.13/
284 & 5X,'scaling
FUNCTION for gas density. . . .=
',I10)
285 1102 FORMAT(
286 & 5X,'iterative newton solver '/)
subroutine hm_get_floatv(name, rval, is_available, lsubmodel, unitab)
subroutine hm_option_is_encrypted(is_encrypted)
for(i8=*sizetab-1;i8 >=0;i8--)
type(alemuscl_param_) alemuscl_param
integer, parameter nchartitle