OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
m11law.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!|| m11law ../engine/source/materials/mat/mat011/m11law.F
25!||--- called by ------------------------------------------------------
26!|| mmain ../engine/source/materials/mat_share/mmain.F90
27!||--- calls -----------------------------------------------------
28!|| finter ../engine/source/tools/curve/finter.F
29!|| m11vs2 ../engine/source/materials/mat/mat011/m11vs2.F
30!|| m11vs3 ../engine/source/materials/mat/mat011/m11vs3.F
31!||--- uses -----------------------------------------------------
32!|| ale_connectivity_mod ../common_source/modules/ale/ale_connectivity_mod.F
33!|| alefvm_mod ../common_source/modules/ale/alefvm_mod.F
34!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
35!||====================================================================
36 SUBROUTINE m11law(PM ,OFF ,SIG ,EINT ,RHO ,
37 2 T ,RK ,RE ,VO ,
38 3 ALE_CONNECT,IX ,IPARG ,ELBUF_TAB,V ,
39 4 X ,EIS ,STIF ,W ,VD2 ,
40 5 VDX ,VDY ,VDZ ,VIS ,MAT ,
41 6 VOL ,QVIS,RHO0 ,BUFVOIS ,IPM ,
42 7 NPF ,TF ,SSP_EQ,MOM ,NEL )
43C-----------------------------------------------
44C M o d u l e s
45C-----------------------------------------------
46 USE elbufdef_mod
47 USE ale_connectivity_mod
48 USE alefvm_mod
49C-----------------------------------------------
50C I m p l i c i t T y p e s
51C-----------------------------------------------
52#include "implicit_f.inc"
53C-----------------------------------------------
54C G l o b a l P a r a m e t e r s
55C-----------------------------------------------
56#include "mvsiz_p.inc"
57C-----------------------------------------------
58C C o m m o n B l o c k s
59C-----------------------------------------------
60#include "com01_c.inc"
61#include "com08_c.inc"
62#include "vect01_c.inc"
63#include "param_c.inc"
64C-----------------------------------------------
65C D u m m y A r g u m e n t s
66C-----------------------------------------------
67 INTEGER IX(*), IPARG(*),MAT(*),IPM(NPROPMI,*),NEL
68 my_real
69 . PM(NPROPM,*), OFF(*), SIG(NEL,6), EINT(*), RHO(*), T(*), RK(*),
70 . RE(*), V(3,*), W(3,*),VIS(*),VOL(*),QVIS(*),
71 . X(3,*), VO(*), EIS(*),STIF(*),VD2(*),VDX(*) ,VDY(*) ,VDZ(*),
72 . RHO0(*),BUFVOIS(*),SSP_EQ(*), MOM(NEL,3)
73 INTEGER :: NPF(*)
74 my_real :: tf(*)
75 my_real ,EXTERNAL :: finter
76 TYPE (ELBUF_STRUCT_), DIMENSION(NGROUP) :: ELBUF_TAB
77 TYPE(t_ale_connectivity), INTENT(IN) :: ALE_CONNECT
78C-----------------------------------------------
79C L o c a l V a r i a b l e s
80C-----------------------------------------------
81 INTEGER I, MX, ITYP, IRFUN, IPFUN, IEFUN, ITFUN, IKFUN,
82 . isfun, inod, ivxfun, ivyfun, ivzfun, ialefvm_flg
83 my_real
84 . gam(mvsiz), p0(mvsiz) , vcrt(mvsiz), gam1(mvsiz), p(mvsiz),
85 . vel(mvsiz), gamrp(mvsiz), rhov(mvsiz), pv(mvsiz) , eiv(mvsiz),
86 . rev(mvsiz), rkv(mvsiz) , tv(mvsiz) , expt(mvsiz),
87 . vxv(mvsiz), vyv(mvsiz) , vzv(mvsiz) ,
88 . fac, cart, e0, xe0, rhof, pf, ef, t0,
89 . xk0, psh, rhon, pn, en, gv, rsr, rv, c1, dc, alp, du, dp,
90 . rhoc,rhoc2,dydx, vxn(mvsiz), vyn(mvsiz), vzn(mvsiz),
91 . vx0(mvsiz), vy0(mvsiz), vz0(mvsiz),
92 . xmomv(mvsiz), ymomv(mvsiz), zmomv(mvsiz) ,
93 . xmomn,ymomn,zmomn
94 my_real
95 . dvp(mvsiz), vn(mvsiz),vimp(mvsiz),ssp,tscal,
96 . rho0_1, gam_1, p0_1, vcrt_1, gamrp_1, gam1_1, vx_1, vy_1, vz_1, mass
97C-----------------------------------------------
98C=======================================================================
99
100 !-----------------------------------!
101 ! READING USER INPUTS !
102 !-----------------------------------!
103 mx = mat(1)
104 rho0_1 = pm( 1,mx)
105 gam_1 = pm(25,mx)
106 p0_1 = pm(31,mx)
107 vcrt_1 = pm(27,mx)
108 gamrp_1 = pm(28,mx)
109 gam1_1 = pm(29,mx)
110 vx_1 = pm(101,mx)
111 vy_1 = pm(102,mx)
112 vz_1 = pm(103,mx)
113 ialefvm_flg = ipm(251,mx)
114 ityp = pm(50,mx)
115 e0 = pm(23,mx)
116 xe0 = pm(33,mx)
117 rhof = pm(35,mx)
118 pf = pm(36,mx)
119 ef = pm(37,mx)
120 irfun = ipm(11,mx)
121 ipfun = ipm(12,mx)
122 iefun = ipm(13,mx)
123 itfun = ipm(14,mx)
124 ikfun = ipm(15,mx)
125 isfun = ipm(16,mx)
126 ivxfun = ipm(18,mx)
127 ivyfun = ipm(19,mx)
128 ivzfun = ipm(20,mx)
129 tscal = pm(40,mx)*tt
130 inod = nint(pm(51,mx))
131 t0 = pm(79,mx)
132 xk0 = pm(87,mx)
133 psh = pm(88,mx)
134
135 !-----------------------------------!
136 ! ELEM INIT. FROM USER INPUTS !
137 !-----------------------------------!
138 DO i=1,nel
139 rho0(i) = rho0_1
140 gam(i) = gam_1
141 p0(i) = p0_1
142 vcrt(i) = vcrt_1
143 gamrp(i) = gamrp_1
144 gam1(i) = gam1_1
145 vx0(i) = vx_1
146 vy0(i) = vy_1
147 vz0(i) = vz_1
148 tv(i) = t0 ! Consequently, if the adjacent material does not have an allocated temperature,
149 ! it will retain its initial value.
150 ENDDO
151
152 DO i=1,nel
153 stif(i) = zero
154 eis(i) = eis(i) + eint(i) / vol(i)
155 ENDDO
156
157 IF(alefvm_param%IEnabled == 1)THEN
158 IF(ityp /= 3)THEN
159 DO i=1,nel
160 IF(ivxfun>0)THEN
161 vxn(i) = vx0(i)*finter(ivxfun,tscal,npf,tf,dydx)
162 ELSE
163 vxn(i) = vx0(i)
164 ENDIF
165 IF(ivyfun>0)THEN
166 vyn(i) = vy0(i)*finter(ivyfun,tscal,npf,tf,dydx)
167 ELSE
168 vyn(i) = vy0(i)
169 ENDIF
170 IF(ivzfun>0)THEN
171 vzn(i) = vz0(i)*finter(ivzfun,tscal,npf,tf,dydx)
172 ELSE
173 vzn(i) = vz0(i)
174 ENDIF
175 ENDDO
176 else!IF(ITYP /= 3)
177 DO i=1,nel
178 mass = rho(i) * vol(i)
179 vxn(i) = mom(i,1) / mass
180 vyn(i) = mom(i,2) / mass
181 vzn(i) = mom(i,3) / mass
182 ENDDO
183 ENDIF
184 ELSE
185 vxn(i) = zero
186 vyn(i) = zero
187 vzn(i) = zero
188 endif!IF(IALEFVM==1)
189
190 !-----------------------------------!
191 ! GET ADJACENT VALUES !
192 !-----------------------------------!
193 IF(n2d==0) THEN
194 CALL m11vs3(pm ,iparg ,ix ,ale_connect ,elbuf_tab ,v ,
195 2 x ,dvp ,vn ,w ,vel ,vd2,
196 3 vdx ,vdy ,vdz ,mat ,rhov ,pv ,
197 4 eiv ,rev ,rkv ,tv ,
198 5 ssp_eq,vxv ,vyv ,vzv ,ialefvm_flg ,
199 6 vxn ,vyn ,vzn ,mom ,rho ,vol,
200 7 xmomv ,ymomv ,zmomv ,bufvois ,nel )
201 fac=eight
202 ELSE
203 CALL m11vs2(pm ,iparg ,ix ,ale_connect ,elbuf_tab ,v ,
204 2 x ,dvp ,vn ,w ,vel ,vd2 ,
205 3 vdy ,vdz ,vis ,mat ,rhov ,pv ,
206 4 eiv ,rev ,rkv ,tv ,
207 5 ssp_eq,vxv ,vyv ,vzv ,ialefvm_flg ,
208 6 vxn ,vyn ,vzn ,mom ,rho ,vol ,
209 7 xmomv ,ymomv ,zmomv ,bufvois ,nel )
210 fac=four
211 ENDIF
212
213 mx = mat(1)
214 DO i=1,nel
215 cart = pm(34,mx)
216 qvis(i) = zero
217 expt(i) = exp(cart*tt**2)
218 ENDDO
219
220
221 !-----------------------------------!
222 ! LOOP ON ELEMENTS !
223 !-----------------------------------!
224 DO i=1,nel
225
226 IF(inod/=0)THEN
227 vel(i) = v(1,inod)**2+v(2,inod)**2+v(3,inod)**2
228 ENDIF
229
230 !---------------------------------!
231 ! IMPOSED STATE ITYP=0,1 !
232 !---------------------------------!
233 IF(ityp<2)THEN
234 rhov(i) = rho0(i)
235 pv(i) = p0(i)
236 eiv(i) = e0
237 vxv(i) = vx0(i)
238 vyv(i) = vy0(i)
239 vzv(i) = vz0(i)
240 !---------------------------------!
241 ! OUTGOING FLUX ITYP=ALL !
242 !---------------------------------!
243 ELSEIF(vn(i)<zero)THEN
244 rhov(i) = rho(i)
245 eiv(i) = eint(i)/vol(i)
246 rkv(i) = rk(i)/vol(i)
247 rev(i) = re(i)/vol(i)
248 tv(i) = t(i)
249 IF(alefvm_param%IEnabled == 1)THEN !otherwise %MOM is not allocated
250 mass = rho(i) * vol(i)
251 vxv(i) = mom(i,1)/mass
252 vyv(i) = mom(i,2)/mass
253 vzv(i) = mom(i,3)/mass
254 xmomv(i)= mom(i,1)
255 ymomv(i)= mom(i,2)
256 zmomv(i)= mom(i,3)
257 ENDIF
258 ENDIF
259
260 IF(irfun>0)THEN
261 rhon = rho0(i)*finter(irfun,tscal,npf,tf,dydx)
262 ELSEIF(irfun<0)THEN
263 rhon = rhof+(rho0(i)-rhof)*expt(i)
264 ELSE
265 rhon = rhov(i)
266 ENDIF
267
268 IF(ipfun>0)THEN
269 pn = p0(i)*finter(ipfun,tscal,npf,tf,dydx)
270 ELSEIF(ipfun<0)THEN
271 pn = pf+(p0(i)-pf)*expt(i)
272 ELSE
273 pn = pv(i)
274 ENDIF
275
276 IF(ityp==5 .OR. ityp==6 )THEN
277 IF(ipfun>0)THEN
278 vimp = p0(i)*finter(ipfun,tscal,npf,tf,dydx)
279 ELSE
280 vimp(i)=p0(i)
281 ENDIF
282 ENDIF
283
284 IF(iefun>0)THEN
285 en = e0*finter(iefun,tscal,npf,tf,dydx)
286 ELSEIF(iefun<0)THEN
287 en =ef+(e0-ef)*expt(i)
288 ELSE
289 en =eiv(i)
290 ENDIF
291
292 IF(jtur/=0)THEN
293 IF(ikfun>0)THEN
294 rk(i) = rho0(i)*finter(ikfun,tscal,npf,tf,dydx)
295 ELSE
296 rk(i) = rkv(i)
297 ENDIF
298 IF(isfun>0)THEN
299 re(i) = rho0(i)*finter(isfun,tscal,npf,tf,dydx)
300 ELSE
301 re(i) = rev(i)
302 ENDIF
303 ENDIF
304
305 IF(itfun>0)THEN
306 t(i) = t0*finter(itfun,tscal,npf,tf,dydx)
307 ELSE
308 t(i) = tv(i)
309 ENDIF
310
311 IF(alefvm_param%IEnabled==1)THEN
312 IF(ivxfun>0)THEN
313 xmomn = vx0(i)*finter(ivxfun,tscal,npf,tf,dydx)*rhon*vol(i)
314 ELSE
315 xmomn = xmomv(i)
316 ENDIF
317 IF(ivyfun>0)THEN
318 ymomn = vy0(i)*finter(ivyfun,tscal,npf,tf,dydx)*rhon*vol(i)
319 ELSE
320 ymomn = ymomv(i)
321 ENDIF
322 IF(ivzfun>0)THEN
323 zmomn = vz0(i)*finter(ivzfun,tscal,npf,tf,dydx)*rhon*vol(i)
324 ELSE
325 zmomn = zmomv(i)
326 ENDIF
327 ENDIF
328
329 !-------------------------------------------!
330 ! CAS 0 et 1 STAGNATION POINT !
331 !-------------------------------------------!
332 !ALEFVM no yet implemented
333 IF(ityp==0)THEN
334C GAZ
335 vel(i) = min(vel(i),vcrt(i))
336 dc=pm(99,mx)
337 gv = dc*gamrp(i)*rhon*vel(i)/(pn+psh)
338 IF(gv>em4)THEN
339 rsr = (one - gv)**gam1(i)
340 p(i) = (pn+psh)*rsr**gam(i)-psh
341 ELSE
342 rv = half*rhon*vel(i)*dc
343 rsr = one - rv/((pn+psh)*gam(i))
344 p(i) = pn - rv
345 ENDIF
346 rho(i) = rsr*rhon
347 p(i) = p(i) + two_third*rk(i)
348 eint(i)= gam1(i)*(p(i)+psh)
349 ssp_eq(i) = sqrt(gam1(i)*(pn+psh)/rhon)
350 ELSEIF(ityp==1)THEN
351C LIQUIDE
352 dc = pm(99,mx)
353 c1 = pm(32,mx)
354 rho(i) = c1*rhon/(c1+half*dc*rhon*vel(i))
355 p(i) = pn-half*dc*rho(i)*vel(i)
356 eint(i) = (one - rho(i)/rhon)*(p(i)+psh)+en
357 ssp_eq(i) = sqrt(c1/rhon)
358 IF(ipfun/=0.AND.jtur/=0)p(i)=p(i)+ two_third*rk(i)
359
360 !-------------------------------------------!
361 ! CAS 2 IMPOSED STATE !
362 !-------------------------------------------!
363 ELSEIF(ityp==2.OR.ityp==8)THEN
364 rho(i) = rhon
365 p(i) = pn
366 eint(i) = en
367 ssp_eq(i) = ep20
368 IF(ipfun/=0.AND.jtur/=0)p(i)=p(i)+two_third*rk(i)
369 IF(alefvm_param%IEnabled==1)THEN
370 mom(i,1) = xmomn
371 mom(i,2) = ymomn
372 mom(i,3) = zmomn
373 ENDIF
374
375 !-------------------------------------------!
376 ! CAS 3 NON REFLECTING BOUNDARY !
377 !-------------------------------------------!
378 ELSEIF(ityp==3 .OR. ityp==7)THEN
379 rho(i) = rhon
380 eint(i) = en
381 rhoc = pm(97,mx)
382 alp = pm(98,mx)*dt1
383 IF(ipfun/=0.AND.jtur/=0)pn=pn + two_third*rk(i)
384 IF(ipfun/=0 .AND. vn(i)<=0)pn=pn - half*rhon*vn(i)**2
385 IF(alefvm_param%IEnabled==1)THEN
386 mom(i,1) = xmomn
387 mom(i,2) = ymomn
388 mom(i,3) = zmomn
389 ENDIF
390C
391 IF(tt>zero)THEN
392 IF(ipfun==0 .OR. vn(i)<=zero .OR. ityp==3)THEN
393 du = (one-alp)*rhoc*(vn(i)-vo(i))
394 ELSE
395 du = (one - alp)*rhoc*((vn(i)-vo(i))-vn(i)*dvp(i))
396 ENDIF
397 ELSE
398 du = zero
399 sig(i,1) = -pn
400 ENDIF
401 dp = alp*(sig(i,1)+pn)
402 vo(i) = vn(i)
403 p(i) = -sig(i,1)+dp+du
404 ssp_eq(i) = rhoc / rho0(i) !adjacent value can also be output but requires SPMD treatment.
405
406 !-------------------------------------------!
407 ! CAS 4 NON REFLECTING BOUNDARY !
408 ! + TRANSVERSE RIGIDITY !
409 !-------------------------------------------!
410 ELSEIF(ityp==4)THEN
411 rho(i) = rhon
412 eint(i)= en
413 rhoc=pm(97,mx)
414 rhoc2=rhon*(rhoc/rho0(i))**2
415 alp=pm(98,mx)*dt1
416 IF(ipfun/=0.AND.jtur/=0)pn=pn + two_third*rk(i)
417 IF(ipfun/=0 .AND. vn(i)<=zero)pn=pn - half*rhon*vn(i)**2
418 IF(tt>zero)THEN
419 du=(one - alp)*(rhoc*(vn(i)-vo(i))-rhoc2*dvp(i))
420 ELSE
421 du=zero
422 sig(i,1)=-pn
423 ENDIF
424
425 dp = alp*(sig(i,1)+pn)
426 vo(i) = vn(i)
427 p(i) = -sig(i,1)+dp+du
428 ssp_eq(i) = rhoc / rho0(i)
429 !-------------------------------------------!
430 ! CAS 5 IMPOSED VELOCITY !
431 ! + NON REFLECTING !
432 !-------------------------------------------!
433 ELSEIF(ityp==5)THEN
434 IF(tt==zero)vo(i) = pv(i)
435 rho(i) = rhon
436 eint(i) = en
437 alp = pm(98,mx)*dt1
438 rhoc = pm(97,mx)
439 vo(i) = vo(i)+alp*(pv(i)-vo(i))
440 p(i) = rhoc*(vn(i)-vimp(i))+vo(i)
441 ssp_eq(i) = rhoc / rho0(i)
442
443 !-------------------------------------------!
444 ! CAS 6 IMPOSED MASS FLOW !
445 ! + NRF !
446 ! + CONSTANT SSP !
447 ! + RHO = RO_ADJACENT !
448 !-------------------------------------------!
449 ELSEIF(ityp==6)THEN
450 IF(tt==zero)vo(i) = pv(i)
451 rho(i) = rhon
452 eint(i) = en
453 alp = pm(98,mx)*dt1
454 ssp = pm(97,mx)/rho0(i)
455 vo(i) = vo(i)+alp*(pv(i)-vo(i))
456 p(i) = ssp*(rhov(i)*vn(i)-vimp(i))+vo(i)
457 ssp_eq(i) = ssp !adjacent value can also be output but requires SPMD treatment.
458 ENDIF
459 ENDDO !next I
460
461C-----------
462 DO i=1,nel
463 eis(i) = eis(i) - eint(i)
464 sig(i,1) = -p(i)
465 sig(i,2) = -p(i)
466 sig(i,3) = -p(i)
467 sig(i,4) = zero
468 sig(i,5) = zero
469 sig(i,6) = zero
470 ENDDO
471C-----------
472 RETURN
473 END
m11law
subroutine m11law(pm, off, sig, eint, rho, t, rk, re, vo, ale_connect, ix, iparg, elbuf_tab, v, x, eis, stif, w, vd2, vdx, vdy, vdz, vis, mat, vol, qvis, rho0, bufvois, ipm, npf, tf, ssp_eq, mom, nel)
Definition m11law.F:43
m11vs2
subroutine m11vs2(pm, iparg, ixq, ale_connect, elbuf_tab, v, x, dvp, vn, w, vel, vd2, vdy, vdz, vis, mat, rhov, pv, eiv, rev, rkv, tv, ssp_eq, vxv, vyv, vzv, ialefvm_flg, vx, vy, vz, mom, rho, vol, xmomv, ymomv, zmomv, bufvois, nel)
Definition m11vs2.F:39
m11vs3
subroutine m11vs3(pm, iparg, ixs, ale_connect, elbuf_tab, v, x, dvp, vn, w, vel, vd2, vdx, vdy, vdz, mat, rhov, pv, eiv, rev, rkv, tv, ssp_eq, vxv, vyv, vzv, ialefvm_flg, vx, vy, vz, mom, rho, vol, xmomv, ymomv, zmomv, bufvois, nel)
Definition m11vs3.F:39
min
#define min(a, b)
Definition macros.h:20
ale_connectivity_mod
Definition ale_connectivity_mod.F:223
alefvm_mod
Definition alefvm_mod.F:88
alefvm_mod::alefvm_param
type(alefvm_param_), target alefvm_param
Definition alefvm_mod.F:121