OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps37.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| sigeps37 ../engine/source/materials/mat/mat037/sigeps37.F
25!||--- called by ------------------------------------------------------
26!|| mulaw ../engine/source/materials/mat_share/mulaw.F90
27!||--- calls -----------------------------------------------------
28!|| finter ../engine/source/tools/curve/finter.F
29!||====================================================================
30 SUBROUTINE sigeps37 (
31 1 NEL ,NUPARAM,NUVAR ,NFUNC ,IFUNC ,NPF ,
32 2 TF ,TIME ,TIMESTEP,UPARAM ,RHO0 ,RHO ,
33 3 VOLUME ,EINT ,
34 4 EPSPXX ,EPSPYY ,EPSPZZ ,EPSPXY ,EPSPYZ ,EPSPZX ,
35 5 DEPSXX ,DEPSYY ,DEPSZZ ,DEPSXY ,DEPSYZ ,DEPSZX ,
36 6 EPSXX ,EPSYY ,EPSZZ ,EPSXY ,EPSYZ ,EPSZX ,
37 7 SIGOXX ,SIGOYY ,SIGOZZ ,SIGOXY ,SIGOYZ ,SIGOZX ,
38 8 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
39 9 SIGVXX ,SIGVYY ,SIGVZZ ,SIGVXY ,SIGVYZ ,SIGVZX ,
40 A SOUNDSP,VISCMAX,UVAR ,OFF ,IX ,NIX ,
41 B NFT )
42C-----------------------------------------------
43C D e s c r i p t i o n
44C-----------------------------------------------
45C Biphasic material law : air / liquid.
46C air is modeled with ideal gas EOS : PV**gamma = constant
47C liquid is modeled with linear EOS : P= P0+C1.mu where C1 is bulk modulus and P0 is initial pressure
48C massic percentage determines the ratio AIR/LIQUID in the cell
49C Purpose is to compute equilibrium Pair=Pwater, this leads to mu_air and mu_water (iterative method, niter=2 is used)
50C
51C ISOLVER = 1 (default) is the legacy solver. Sound Speed is computed from water whatever is the mixture. There is only 2 iterations.
52C ISOLVER = 2 is an alternative solver using a Newton Algorithm with (max 20 iterations). Convergence criteria introduced and sound speed computed from mixture.
53C This sound speed formulation may leads to global value potentially lower than each submaterial value.
54C recommended : Courant number <= 0.5 (dt_scale)
55C
56C Output
57C !UVAR(I,1) : massic percentage of liquid * global density (rho1*V1/V : it needs to give liquid mass multiplying by element volume in aleconve.F)
58C !UVAR(I,2) : density of gas
59C !UVAR(I,3) : density of liquid
60C !UVAR(I,4) : volumetric fraction of liquid
61C !UVAR(I,5) : volumetric fraction of gas
62C
63C !---------+---------+---+---+--------------------------------------------
64C ! VAR | SIZE |TYP| RW| DEFINITION
65C !---------+---------+---+---+--------------------------------------------
66C ! NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
67C ! NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
68C ! NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
69C !---------+---------+---+---+--------------------------------------------
70C ! NFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW
71C ! IFUNC | NFUNC | I | R | FUNCTION INDEX
72C ! NPF | * | I | R | FUNCTION ARRAY
73C ! TF | * | F | R | FUNCTION ARRAY
74C !---------+---------+---+---+--------------------------------------------
75C ! TIME | 1 | F | R | CURRENT TIME
76C ! TIMESTEP| 1 | F | R | CURRENT TIME STEP
77C ! UVAR | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
78C ! RHO0 | NEL | F | R | INITIAL DENSITY
79C ! RHO | NEL | F | R | DENSITY
80C ! VOLUME | NEL | F | R | VOLUME
81C ! EINT | NEL | F | R | TOTAL INTERNAL ENERGY
82C ! EPSPXX | NEL | F | R | STRAIN RATE XX
83C ! EPSPYY | NEL | F | R | STRAIN RATE YY
84C ! ... | | | |
85C ! DEPSXX | NEL | F | R | STRAIN INCREMENT XX
86C ! DEPSYY | NEL | F | R | STRAIN INCREMENT YY
87C ! ... | | | |
88C ! EPSXX | NEL | F | R | STRAIN XX
89C ! EPSYY | NEL | F | R | STRAIN YY
90C ! ... | | | |
91C ! SIGOXX | NEL | F | R | OLD ELASTO PLASTIC STRESS XX
92C ! SIGOYY | NEL | F | R | OLD ELASTO PLASTIC STRESS YY
93C ! ... | | | |
94C !---------+---------+---+---+--------------------------------------------
95C ! SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
96C ! SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
97C ! ... | | | |
98C ! SIGVXX | NEL | F | W | VISCOUS STRESS XX
99C ! SIGVYY | NEL | F | W | VISCOUS STRESS YY
100C ! ... | | | |
101C ! SOUNDSP | NEL | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
102C ! VISCMAX | NEL | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
103C !---------+---------+---+---+--------------------------------------------
104C ! UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
105C ! OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
106C !---------+---------+---+---+--------------------------------------------
107C
108C-----------------------------------------------
109C I m p l i c i t T y p e s
110C-----------------------------------------------
111#include "implicit_f.inc"
112C-----------------------------------------------
113C I N P U T A r g u m e n t s
114C-----------------------------------------------
115 INTEGER NEL, NUPARAM, NUVAR
116 my_real TIME,TIMESTEP,UPARAM(NUPARAM),
117 . RHO(NEL),RHO0(NEL),VOLUME(NEL),EINT(NEL),
118 . EPSPXX(NEL),EPSPYY(NEL),EPSPZZ(NEL),
119 . EPSPXY(NEL),EPSPYZ(NEL),EPSPZX(NEL),
120 . DEPSXX(NEL),DEPSYY(NEL),DEPSZZ(NEL),
121 . DEPSXY(NEL),DEPSYZ(NEL),DEPSZX(NEL),
122 . EPSXX(NEL) ,EPSYY(NEL) ,EPSZZ(NEL) ,
123 . EPSXY(NEL) ,EPSYZ(NEL) ,EPSZX(NEL) ,
124 . SIGOXX(NEL),SIGOYY(NEL),SIGOZZ(NEL),
125 . SIGOXY(NEL),SIGOYZ(NEL),SIGOZX(NEL)
126C-----------------------------------------------
127C O U T P U T A r g u m e n t s
128C-----------------------------------------------
129 my_real, INTENT(INOUT) ::
130 . signxx(nel),signyy(nel),signzz(nel),
131 . signxy(nel),signyz(nel),signzx(nel),
132 . sigvxx(nel),sigvyy(nel),sigvzz(nel),
133 . sigvxy(nel),sigvyz(nel),sigvzx(nel),
134 . soundsp(nel),viscmax(nel)
135C-----------------------------------------------
136C I N P U T O U T P U T A r g u m e n t s
137C-----------------------------------------------
138 my_real,INTENT(INOUT) :: uvar(nel,nuvar), off(nel)
139 INTEGER, INTENT(IN) :: NPF(*), NFUNC, IFUNC(NFUNC), NIX, IX(NIX,*), NFT
140 my_real, INTENT(IN) :: tf(*)
141 my_real, EXTERNAL :: finter
142C-----------------------------------------------
143C L o c a l V a r i a b l e s
144C-----------------------------------------------
145 INTEGER I,J, ISOLVER, NITER, ITER
146 my_real
147 . SSP,VIS,VIS2,VIS3,VV,C1,C2,C12,R1,R2,PMIN,A2,RHO10,RHO20,
148 . RHO1,RHO2,A1,A,B,C, B1,B2,P,GAM,P0,GPR,POLD,
149 . pn2,rhn2,visa1,visb1,visa2,visb2,dydx,rhoscale,tol, vol, mas, mas1, mas2,
150 . rho1t,rho2t, error, p1,p2,f1,f2,df11,df12,df21,df22,det,drho1,drho2,
151 . mu1p1, mu2p1, psh, ssp1, ssp2
152C-----------------------------------------------
153C S o u r c e C o d e
154C-----------------------------------------------
155
156 !------------------------------------!
157 ! USER PARAMETERS !
158 !------------------------------------!
159 visa1 = uparam(1)
160 visb1 = uparam(3)
161 visa2 = uparam(13)
162 visb2 = uparam(15)
163 c1 = uparam(4)
164 gam = uparam(5)
165 r1 = uparam(6)
166 gpr = uparam(7)
167 pmin = uparam(8)
168 p0 = uparam(9)
169 rho10 = uparam(11)
170 rho20 = uparam(12)
171 psh = uparam(16)
172 isolver = nint(uparam(17))
173 rhoscale = one
174 IF(ifunc(1)>0)rhoscale=finter(ifunc(1),time,npf,tf,dydx)
175
176 !------------------------------------!
177 ! INITIALIZATION TIME==ZERO !
178 !------------------------------------!
179 IF(time==zero)THEN
180 DO i=1,nel
181 p = max(em30,(-sigoxx(i)-sigoyy(i)-sigozz(i))*third)-psh
182 IF(gam*c1>=em30)THEN !if Liquid and gas correctly defined
183 mu1p1 = ((p-p0)/c1+one)
184 mu2p1 =( p/p0)**(one/gam)
185 rho1 = rho10*mu1p1
186 rho2 = rho20*mu2p1
187 a = (rho(i)-rho2)/(rho1-rho2)
188 uvar(i,1) = a*rho1
189 uvar(i,2) = rho2
190 uvar(i,3) = rho1
191 uvar(i,4) = a
192 IF(uvar(i,4)<em20)uvar(i,4)=zero
193 uvar(i,5) = one-uvar(i,4)
194 ELSE !boundary element
195 uvar(i,3)=rho(i)
196 ENDIF
197 ENDDO
198 ELSE
199 DO i=1,nel
200 uvar(i,1) = uvar(i,1)/volume(i) !updated in aleconve.F
201 ENDDO
202 ENDIF
203
204 !------------------------------------!
205 ! CASE OF BOUNDARY ELEMENT INPUT !
206 ! C1=0:not a fluid ; gam=0:not a gas !
207 !------------------------------------!
208 IF(gam*c1<em30)THEN
209 DO i=1,nel
210 IF(uvar(i,3)/rho10<half)THEN
211 uvar(i,1) = zero !massic percentage of liquid
212 rho(i) = uvar(i,3)*rhoscale
213 uvar(i,2) = rho(i)
214 uvar(i,4) = zero ! no mass of liquid then no volume
215 uvar(i,5) = one
216 ELSE
217 uvar(i,1) = uvar(i,3)
218 uvar(i,2) = rho20
219 rho(i) = uvar(i,3)
220 uvar(i,4) = one
221 uvar(i,5) = zero
222 ENDIF
223 soundsp(i) = em30
224 signxx(i) = sigoxx(i)
225 signyy(i) = sigoyy(i)
226 signzz(i) = sigozz(i)
227 signxy(i) = sigoxy(i)
228 signyz(i) = sigoyz(i)
229 signzx(i) = sigozx(i)
230 ENDDO
231 RETURN
232 ENDIF
233
234 IF(isolver==2)THEN
235 !------------------------------------!
236 ! ALTERNATIVE SOLVER (Newton) !
237 !------------------------------------!
238 tol = em10
239 niter = 20
240 DO i=1,nel
241 vol = volume(i)
242 mas = rho(i) * vol
243 mas1 = uvar(i, 1) * vol
244 mas2 = mas - mas1
245 rho2 = uvar(i, 2)
246 rho1 = uvar(i, 3)
247 rho1t = rho1 / rho10
248 rho2t = rho2 / rho20
249 pold = p0 * rho2t**gam
250 IF (mas1 / mas < em10) THEN
251 !! Phase 2 only is present
252 uvar(i, 1) = zero
253 uvar(i, 4) = zero
254 uvar(i, 5) = one
255 rho2 = mas / vol
256 uvar(i, 2) = rho2
257 p = p0 * (rho2/rho20)**gam
258 ELSEIF (mas2 / mas < em10) THEN
259 !! Phase 1 only is present
260 rho1 = mas / vol
261 uvar(i, 1) = rho1
262 uvar(i, 3) = rho1
263 uvar(i, 4) = one
264 uvar(i, 5) = zero
265 p = r1 * rho1 - c1 + p0
266 ELSE
267 error = ep30
268 iter = 1
269 DO WHILE(iter < niter .AND. error > tol)
270 p1 = r1 * rho1 - c1 + p0
271 p2 = p0 * (rho2/rho20)**gam
272 f1 = mas1 / rho1 + mas2 / rho2 - vol
273 f2 = p1 - p2
274 df11 = - mas1 / (rho1 * rho1)
275 df12 = - mas2 / (rho2 * rho2)
276 df21 = r1
277 df22 = - gam * p0 / (rho20**gam) * rho2**(gam - one)
278 det = df11 * df22 - df12 * df21
279 drho1 = (-df22 * f1 + df12 * f2) / det
280 drho2 = (df21 * f1 - df11 * f2) / det
281 drho1 = min(three * rho1, max(drho1, - half * rho1))
282 drho2 = min(three * rho2, max(drho2, - half * rho2))
283 rho1 = rho1 + drho1
284 rho2 = rho2 + drho2
285 error = abs(drho1 / rho1) + abs(drho2 / rho2)
286 iter = iter + 1
287 ENDDO
288 IF (error > tol) THEN
289 print*, "*** WARNING LAW37, convergence tolerance ", error, tol
290 ENDIF
291 p = r1 * rho1 - c1 + p0
292 ENDIF
293 ssp1 = r1 * rho1
294 ssp2 = gam * p0 * (rho2/rho20)**gam
295 b1 = uvar(i, 1)
296 b2 = rho(i) - b1
297 !---output
298 uvar(i,2) = rho2
299 uvar(i,3) = rho1
300 uvar(i,4) = uvar(i,1)/rho1
301 IF(uvar(i,4)<em20)uvar(i,4)=zero
302 uvar(i,5) = one-uvar(i,4)
303 IF (ssp1 > zero) THEN
304 ssp1 = uvar(i,4) / ssp1
305 ELSE
306 ssp1 = zero
307 ENDIF
308 IF (ssp2 > zero) THEN
309 ssp2 = uvar(i,5) / ssp2
310 ELSE
311 ssp2 = zero
312 ENDIF
313 ssp = ssp1 + ssp2
314 ssp = sqrt(one / ssp / rho(i))
315 p = max(pmin, p) + psh
316 signxx(i) = -p
317 signyy(i) = -p
318 signzz(i) = -p
319 vis = (b1*rho1*visa1 + b2*rho2*visa2)/rho(i)
320 vis2 = two*vis
321 vis3 = (b1*rho1*visb1 + b2*rho2*visb2)/rho(i)
322 vv = vis3*(epspxx(i)+epspyy(i)+epspzz(i))
323 sigvxx(i) = vis2*epspxx(i)+vv
324 sigvyy(i) = vis2*epspyy(i)+vv
325 sigvzz(i) = vis2*epspzz(i)+vv
326 sigvxy(i) = vis *epspxy(i)
327 sigvyz(i) = vis *epspyz(i)
328 sigvzx(i) = vis *epspzx(i)
329 soundsp(i) = ssp
330 viscmax(i) = vis2 + vis3
331 uvar(i,1) = max(zero, uvar(i,1))
332 uvar(i,2) = max(zero, uvar(i,2))
333 uvar(i,3) = max(zero, uvar(i,3))
334 uvar(i,4) = max(zero, uvar(i,4))
335 uvar(i,5) = max(zero, uvar(i,5))
336 enddo!next I
337
338 ELSE !ISOLVER
339 !------------------------------------!
340 ! LEGACY SOLVER (default) !
341 !------------------------------------!
342 DO i=1,nel
343 rho2 = uvar(i,2)
344 !---iter 1
345 pold = p0 * (rho2/rho20)**gam
346 r2 = gam * pold / rho2
347 c2 = - (one-gam)*pold + p0
348 c12 = c1 - c2
349 b1 = uvar(i,1)
350 b2 = rho(i) - b1
351 a = r1
352 b = half*(b1*r1+b2*r2+c12)
353 c = b1*c12
354 rho1 = ( b + sqrt(max(zero,b*b - a*c)) ) / a
355 p = r1*rho1 - c1
356 rhn2 = max(em30,(p + c2) / r2)
357 !---iter 2
358 pn2 = (pold + p0 * (rhn2/rho20)**gam)
359 r2 = gam * pn2 / (rho2+rhn2)
360 b = half*(b1*r1+b2*r2+c12)
361 rho1 = ( b + sqrt(max(zero,b*b - a*c)) ) / a
362 p = r1*rho1 - c1
363 rho2 = max(em30,(p + c2) / r2)
364 !---output
365 uvar(i,2) = rho2
366 uvar(i,3) = rho1
367 uvar(i,4) = uvar(i,1)/rho1
368 IF(uvar(i,4)<em20)uvar(i,4)=zero
369 uvar(i,5) = one-uvar(i,4)
370 p = max(pmin, p) +p0+psh !P is here a relative pressure
371 signxx(i) = -p
372 signyy(i) = -p
373 signzz(i) = -p
374 vis = (b1*rho1*visa1 + b2*rho2*visa2)/rho(i)
375 vis2 = two*vis
376 vis3 = (b1*rho1*visb1 + b2*rho2*visb2)/rho(i)
377 vv = vis3*(epspxx(i)+epspyy(i)+epspzz(i))
378 sigvxx(i) = vis2*epspxx(i)+vv
379 sigvyy(i) = vis2*epspyy(i)+vv
380 sigvzz(i) = vis2*epspzz(i)+vv
381 sigvxy(i) = vis *epspxy(i)
382 sigvyz(i) = vis *epspyz(i)
383 sigvzx(i) = vis *epspzx(i)
384 soundsp(i) = sqrt(c1/rho1)
385 viscmax(i) = vis2 + vis3
386 uvar(i,1) = max(zero, uvar(i,1))
387 uvar(i,2) = max(zero, uvar(i,2))
388 uvar(i,3) = max(zero, uvar(i,3))
389 uvar(i,4) = max(zero, uvar(i,4))
390 uvar(i,5) = max(zero, uvar(i,5))
391 enddo!next I
392 ENDIF
393
394C-----------------------------------------------
395 RETURN
396 END
sigeps37
subroutine sigeps37(nel, nuparam, nuvar, nfunc, ifunc, npf, tf, time, timestep, uparam, rho0, rho, volume, eint, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, depsxx, depsyy, depszz, depsxy, depsyz, depszx, epsxx, epsyy, epszz, epsxy, epsyz, epszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, sigvxx, sigvyy, sigvzz, sigvxy, sigvyz, sigvzx, soundsp, viscmax, uvar, off, ix, nix, nft)
Definition sigeps37.F:42
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21