OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mqviscb.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!|| mqviscb ../engine/source/materials/mat_share/mqviscb.F
25!||--- called by ------------------------------------------------------
26!|| m1law ../engine/source/materials/mat/mat001/m1law.F
27!|| m1lawi ../engine/source/materials/mat/mat001/m1lawi.F
28!|| m1lawtot ../engine/source/materials/mat/mat001/m1lawtot.F
29!|| m22law ../engine/source/materials/mat/mat022/m22law.F
30!|| m24law ../engine/source/materials/mat/mat024/m24law.F
31!|| m2law ../engine/source/materials/mat/mat002/m2law.F
32!|| m46law ../engine/source/materials/mat/mat046/m46law.f
33!|| mmain ../engine/source/materials/mat_share/mmain.f90
34!|| mulaw ../engine/source/materials/mat_share/mulaw.F90
35!|| sboltlaw ../engine/source/elements/solid/solide/sboltlaw.F
36!|| usermat_solid ../engine/source/materials/mat_share/usermat_solid.F
37!||--- uses -----------------------------------------------------
38!|| ale_mod ../common_source/modules/ale/ale_mod.F
39!|| glob_therm_mod ../common_source/modules/mat_elem/glob_therm_mod.F90
40!|| i22bufbric_mod ../common_source/modules/interfaces/cut-cell-search_mod.F
41!|| i22tri_mod ../common_source/modules/interfaces/cut-cell-search_mod.F
42!||====================================================================
43 SUBROUTINE mqviscb(
44 1 PM, OFF, RHO, RK,
45 2 TEMP, SSP, RE, STI,
46 3 DT2T, NELTST, ITYPTST, AIRE,
47 4 OFFG, GEO, PID, VOL,
48 5 VD2, DELTAX, VIS, D1,
49 6 D2, D3, PNEW, PSH,
50 7 MAT, NGL, QVIS, SSP_EQ,
51 8 VOL0, MSSA, DMELS, IGEO,
52 9 FACQ0, CONDE, DTEL, G_DT,
53 A IPM, RHOREF, RHOSP, NEL,
54 B ITY, ISMSTR, JTUR, JTHE,
55 C JSMS, NPG , GLOB_THERM)
56C============================================================================
57C-----------------------------------------------
58C M o d u l e s
59C-----------------------------------------------
60 USE i22tri_mod
61 USE i22bufbric_mod
62 USE ale_mod
63 use glob_therm_mod
64C-----------------------------------------------
65C I m p l i c i t T y p e s
66C-----------------------------------------------
67#include "implicit_f.inc"
68#include "comlock.inc"
69C-----------------------------------------------
70C G l o b a l P a r a m e t e r s
71C-----------------------------------------------
72#include "mvsiz_p.inc"
73C-----------------------------------------------
74C C o m m o n B l o c k s
75C-----------------------------------------------
76#include "com01_c.inc"
77#include "com04_c.inc"
78#include "com08_c.inc"
79#include "cong1_c.inc"
80#include "impl1_c.inc"
81#include "param_c.inc"
82#include "scr02_c.inc"
83#include "scr07_c.inc"
84#include "scr17_c.inc"
85#include "scr18_c.inc"
86#include "sms_c.inc"
87#include "units_c.inc"
88#include "inter22.inc"
89C-----------------------------------------------
90C D u m m y A r g u m e n t s
91C-----------------------------------------------
92 INTEGER, INTENT(IN) :: NEL
93 INTEGER, INTENT(IN) :: ITY
94 INTEGER, INTENT(IN) :: ISMSTR
95 INTEGER, INTENT(IN) :: JTUR
96 INTEGER, INTENT(IN) :: JTHE
97 INTEGER, INTENT(IN) :: JSMS,NPG
98 INTEGER NELTST,ITYPTST,PID(*),MAT(*), NGL(*),IGEO(NPROPGI,NUMGEO),G_DT,IPM(NPROPMI,NUMMAT)
99 my_real DT2T
100 my_real PM(NPROPM,NUMMAT), OFF(*), RHO(*), RK(*), TEMP(*), RE(*),STI(*),
101 . OFFG(*),GEO(NPROPG,NUMGEO),
102 . VOL(*), VD2(*), DELTAX(*), SSP(*), AIRE(*), VIS(*),
103 . PSH(*), PNEW(*),QVIS(*) ,SSP_EQ(*),CONDE(*),
104 . d1(*), d2(*), d3(*), vol0(*), mssa(*), dmels(*),facq0
105 my_real, INTENT(IN) :: rhoref(mvsiz), rhosp(mvsiz)
106 my_real, INTENT(INOUT) :: dtel(mvsiz)
107 type (glob_therm_) ,intent(inout) :: glob_therm
108C-----------------------------------------------
109C L o c a l V a r i a b l e s
110C-----------------------------------------------
111 INTEGER I,J, MT, K,ICOUNT,LIST(MVSIZ),ERROR, ALE_OR_EULER, ID_MIN
112 my_real dd(mvsiz), al(mvsiz),
113 . dtx(mvsiz), ad(mvsiz), qx(mvsiz), cx(mvsiz), rho0(mvsiz), nrho(mvsiz),
114 . dty(mvsiz), qa, qb, visi, facq,qaa,
115 . cns1, cns2, sph, ak1, bk1, ak2, bk2, tli, akk, xmu, tmu, rpr,
116 . atu, qad, qbd, qaap, dt, nui
117 my_real tidt,tvol,trho,taire, dtmin(mvsiz), rhos(mvsiz),facpg
118 my_real :: cns1_0,cns2_0,qaa_0
119#ifdef MYREAL8
120 my_real, PARAMETER :: real_three = 3.0d0
121 my_real, PARAMETER :: real_one = 1.0d0
122#else
123 my_real, PARAMETER :: real_three = 3.0
124 my_real, PARAMETER :: real_one = 1.0
125#endif
126
127C=======================================================================
128
129 ! by default don't apply /DT/NODA/* to ALE/EULER cells ; unless if /DT/NODA/ALEON is enabled (hidden card introduced for backward compatibility)
130 mt = mat(1)
131 ale_or_euler = 0
132 IF(nint(pm(72,mt))==1 .OR. nint(pm(72,mt))==2)ale_or_euler = 1
133 IF( ale%GLOBAL%I_DT_NODA_ALE_ON==1) ale_or_euler = 0
134
135 IF(impl==zero)THEN
136 DO i=1,nel
137 dd(i)=-d1(i)-d2(i)-d3(i)
138 ad(i)=zero
139 al(i)=zero
140 cx(i)=ssp(i) + sqrt(vd2(i))
141 ENDDO
142 IF((impl_s>0.AND.idyna==0).OR.ismdisp>0)THEN
143 visi=zero
144 facq=zero
145 ELSE
146 visi=one
147 facq=facq0
148 ENDIF
149 IF(imconv<0) THEN
150 DO i=1,nel
151 vol(i)= abs(vol(i))
152 ENDDO
153 ENDIF
154 ELSE
155 DO i=1,nel
156 dd(i)=-d1(i)-d2(i)-d3(i)
157 ad(i)=zero
158 al(i)=zero
159 cx(i)=sqrt(vd2(i))
160 ENDDO
161 visi=zero
162 facq=zero
163 ENDIF
164C
165 IF(n2d>0) THEN
166 facpg=one
167 DO i=1,nel
168 IF(off(i)==one)THEN
169 al(i)=sqrt(aire(i))
170 ad(i)= max(zero,dd(i))
171 ENDIF
172 ENDDO
173 ELSE
174 facpg = npg
175 IF (facpg>one) facpg=facpg**(real_one/real_three)
176 DO i=1,nel
177 IF(off(i)==one)THEN
178 IF(vol(i) > zero)THEN
179 al(i) = vol(i)**(real_one/real_three)
180 ELSE
181 al(i)=zero
182 END IF
183 ad(i)= max(zero,dd(i))
184 ENDIF
185 ENDDO
186 ENDIF
187C
188 !Numerical viscosity is managed inside sigeps51.f for law51
189 !qa, qb should be 0 here otherwise viscosity will be twice.
190 IF(.NOT.(idtmins==2.AND.jsms/=0))THEN
191 mt = mat(1)
192 DO i=1,nel
193 rho0(i)=pm(1,mt)
194 nrho(i)=sqrt(rhoref(i)*rho0(i))
195 ENDDO
196 qa =facq*geo(14,pid(1))
197 qb =facq*geo(15,pid(1))
198 cns1_0=facpg*geo(16,pid(1))
199 cns2_0=facpg*geo(17,pid(1))
200 psh(1:nel)=pm(88,mt)
201 pnew(1:nel)=zero
202 qaa_0 = qa*qa
203 dtmin(1:nel) = geo(172,pid(1))
204 DO i=1,nel
205 cns1=cns1_0*al(i)*nrho(i)*ssp(i)*off(i)
206 cns2=cns2_0*al(i)*nrho(i)*ssp(i)*off(i)
207C
208C warns : Stability vs artificial viscosity should use
209C RHO(I)*AL(I)*(QAA*AL(I)+QB*SSP(I)) / MAX(EM20,RHOREF(I)*DELTAX(I))
210C Same for VIS(I) since SOUNDSPEED is generally computed wrt RHO0, not RHO ../..
211 qaa = qaa_0*ad(i)
212 qx(i)=qb*ssp(i)+al(i) * qaa
213 . + visi*two*vis(i) / max(em20,rho(i)*deltax(i))
214 . + (cns1+visi*cns2) / max(em20,rhoref(i)*deltax(i))
215 qvis(i)=rho(i)*ad(i)*al(i)*(qaa*al(i)+qb*ssp(i))
216 ENDDO
217 ELSE
218 mt = mat(1)
219 DO i=1,nel
220 rho0(i)=pm(1,mt)
221 nrho(i)=sqrt(rhoref(i)*rho0(i))
222 ENDDO
223
224 qa =facq*geo(14,pid(1))
225 qb =facq*geo(15,pid(1))
226 psh(1:nel)=pm(88,mt)
227 pnew(1:nel)=zero
228 qaa_0 = qa*qa
229 dtmin(1:nel) = geo(172,pid(1))
230
231 IF(geo(16,pid(1)) >= zero)THEN
232 cns1_0=facpg*geo(16,pid(1))
233 cns2_0=facpg*geo(17,pid(1))
234 DO i=1,nel
235 cns1=cns1_0*al(i)*nrho(i)*ssp(i)*off(i)
236 cns2=cns2_0*al(i)*nrho(i)*ssp(i)*off(i)
237 qaa = qaa_0*ad(i)
238 nui =rho(i)*al(i)*(qaa*al(i)+qb*ssp(i))
239 qvis(i)=nui*ad(i)
240 qx(i) =(nui + cns1+visi*(two*vis(i)+cns2))/max(em20,rhoref(i)*deltax(i))
241 ENDDO
242 ELSE
243 cns1_0=abs(geo(16,pid(1)))
244 cns2_0=abs(geo(17,pid(1)))
245 DO i=1,nel
246 cns1=cns1_0*nrho(i)*ssp(i)**2*off(i)
247 cns2=cns2_0*nrho(i)*ssp(i)**2*off(i)
248 qaa = qaa_0*ad(i)
249 nui =rho(i)*al(i)*(qaa*al(i)+qb*ssp(i))
250 qvis(i)=nui*ad(i)
251 qx(i) =(nui + cns1+visi*(two*vis(i)+cns2))/max(em20,rhoref(i)*deltax(i))
252 ENDDO
253 END IF
254 ENDIF
255C
256 DO i=1,nel
257 ssp_eq(i) = max(em20,qx(i)+sqrt(qx(i)*qx(i)+cx(i)*cx(i)))
258 dtx(i) = deltax(i) / ssp_eq(i)
259 ENDDO
260C
261 !This would output CFL value without scale factor.
262 !IF (G_DT>0)THEN
263 ! DO I=1,NEL
264 ! DTEL(I) = DTX(I)
265 ! ENDDO
266 !ENDIF
267C
268 IF(int22 > 0)THEN
269 !check if at least one brick has time step smaller than cinematic one => elementary time step.
270 DO i=1,nel
271 IF(dtx(i)<dt22_min)THEN
272 idt_int22 = 0
273 !NELTST = NGL(I)
274 !ITYPTST = ITY
275 EXIT ! just set IDT_INT22
276 ENDIF
277 ENDDO
278 ENDIF
279C
280 IF(jthe > 0)THEN
281 mt = mat(1)
282 sph = pm(69,mt)
283 ak1 = pm(75,mt)
284 bk1 = pm(76,mt)
285 ak2 = pm(77,mt)
286 bk2 = pm(78,mt)
287 tli = pm(80,mt)
288 DO i=1,nel
289 IF(temp(i)<tli)THEN
290 akk=ak1+bk1*temp(i)
291 ELSE
292 akk=ak2+bk2*temp(i)
293 ENDIF
294 IF(jtur/=0)THEN
295 xmu = rho(i)*pm(24,mt)
296 tmu = pm(81,mt)
297 rpr = pm(95,mt)
298 atu=rpr*tmu*rk(i)*rk(i)/(max(em15,re(i)*vol(i))*xmu)
299 akk=akk*(1.+atu)
300 ENDIF
301 dtx(i) = min(dtx(i),half*deltax(i)*deltax(i)*sph/max(akk,em20))
302 ENDDO
303 ENDIF
304C
305 IF(.NOT.(idtmins==2.AND.jsms/=0))THEN
306 IF(kdtsmstr==1.AND.ismstr==1.OR.
307 . ((ismstr==2.OR.ismstr==12).AND.idtmin(1)==3))THEN
308 ! KDTSMSTR==1 en version 5, par defaut.
309 mt = mat(1)
310 DO i=1,nel
311 rho0(i)=pm(1,mt)
312 END DO
313 DO 50 i=1,nel
314 sti(i) = zero
315 IF(off(i)==zero.OR.offg(i)<zero) GO TO 50
316 IF(ipm(251,mt)/=0) GO TO 50
317 IF(n2d==0) THEN
318 tidt = one/dtx(i)
319 IF(offg(i)>one)THEN
320 trho = rho0(i) * tidt
321 tvol = vol0(i) * tidt
322 ELSE
323 trho = rho(i) * tidt
324 tvol = vol(i) * tidt
325 END IF
326 ! STI will be changed to 2*STI/NNE in SxCUMU
327 ! [ STI(I) = FOURTH * TRHO * TVOL for 8node bricks ]
328 sti(i) = trho * tvol
329 ELSE
330 ! small strain is not available in 2D analysis
331 tidt = one/dtx(i)
332 trho = rho(i) * tidt
333 taire = aire(i) * tidt
334 sti(i) = half * trho * taire
335 ENDIF
336C dt2 remplace par dt2t
337 50 CONTINUE
338 IF(ale_or_euler==0)THEN
339 DO i=1,nel
340 dtx(i)= dtfac1(ity)*dtx(i)
341 ENDDO
342 ELSE
343 DO i=1,nel
344 dtx(i)= dtfac1(102)*dtx(i)
345 ENDDO
346 ENDIF
347 !/DT/NODA & /DT/NODA/CST HAS NO EFFECT WITH ALE/EULER
348 IF(ale_or_euler==1 .OR. nodadt==0)THEN
349 DO i=1,nel
350 IF(vol(i)>zero .AND. (off(i)/=zero.AND.offg(i)>=zero))dt2t= min(dtx(i),dt2t)
351 ENDDO
352 ENDIF
353
354 ELSE
355 IF(ismstr == 11) THEN
356 IF(n2d == 0) THEN
357 DO i=1,nel
358 sti(i) = zero
359 IF(off(i)==zero.OR.offg(i)<zero) cycle
360 IF(ipm(251,mt)/=0) cycle
361 tidt = one/dtx(i)
362 trho = rho0(i) * tidt
363 tvol = vol0(i) * tidt
364 sti(i) = trho * tvol
365 ENDDO
366 ELSE
367 DO i=1,nel
368 sti(i) = zero
369 IF(off(i)==zero.OR.offg(i)<zero) cycle
370 IF(ipm(251,mt)/=0) cycle
371 tidt = one/dtx(i)
372 trho = rho(i) * tidt
373 taire = aire(i) * tidt
374 sti(i) = half * trho * taire
375 ENDDO
376 ENDIF ! N2D
377 ELSE
378 DO 60 i=1,nel
379 sti(i) = zero
380 IF(off(i)==zero.OR.offg(i)<zero) GO TO 60
381 IF(ipm(251,mt)/=0) GO TO 60
382 IF(n2d==0) THEN
383 tidt = one/dtx(i)
384 trho = rho(i) * tidt
385 tvol = vol(i) * tidt
386 sti(i) = trho * tvol
387C STI will be changed to 2*STI/NNE in SxCUMU
388C [ STI(I) = FOURTH * TRHO * TVOL for 8node bricks ]
389 ELSE
390 tidt = one/dtx(i)
391 trho = rho(i) * tidt
392 taire = aire(i) * tidt
393 sti(i) = half * trho * taire
394 ENDIF
395C dt2 remplace par dt2t
396 60 CONTINUE
397 ENDIF ! ISMSTR == 11
398 IF(ale_or_euler==0)THEN
399 DO i=1,nel
400 dtx(i)= dtfac1(ity)*dtx(i)
401 ENDDO
402 ELSE
403 DO i=1,nel
404 dtx(i)= dtfac1(102)*dtx(i)
405 ENDDO
406 ENDIF
407 !/DT/NODA & /DT/NODA/CST HAS NO EFFECT WITH ALE/EULER
408 IF(ale_or_euler==1 .OR. nodadt==0)THEN
409 DO i=1,nel
410 IF(vol(i)>zero .AND. (off(i)/=zero.AND.offg(i)>=zero))dt2t= min(dtx(i),dt2t)
411 ENDDO
412 ENDIF
413 END IF
414C
415 ELSE ! IDTMINS == 2 .AND. JSMS /= 0
416 DO i=1,nel
417 dty(i)= dtx(i)
418 dtx(i)= dtfac1(ity)*dtx(i)
419 END DO
420 END IF
421C----
422 IF(imconv==1)THEN
423 IF(idtmin(ity)==1)THEN
424 error=0
425 DO 70 i=1,nel
426 IF(dtx(i)>dtmin1(ity).OR.off(i)==zero.
427 . or.offg(i)<zero) GO TO 70
428 error=1
429 70 CONTINUE
430
431 IF (error==1) THEN
432 tstop = tt
433 DO i=1,nel
434 IF(dtx(i)>dtmin1(ity).OR.off(i)==zero.
435 . or.offg(i)<zero) cycle
436#include "lockon.inc"
437 WRITE(iout,*)
438 . ' **ERROR : TIME STEP LESS OR EQUAL DTMIN FOR SOLID ELEMENT',
439 . ngl(i)
440#include "lockoff.inc"
441 END DO
442#include "lockon.inc"
443 WRITE(istdo,*)
444 . ' **ERROR : TIME STEP LESS OR EQUAL DTMIN FOR SOLID ELEMENTS'
445#include "lockoff.inc"
446 ENDIF
447 ELSEIF(idtmin(ity)==2)THEN
448 icount=0
449 DO 75 i=1,nel
450 IF(dtmin(i)/=zero) THEN
451 IF(dtx(i)>dtmin(i).OR.off(i)==zero.
452 . or.offg(i)<zero) GO TO 75
453 ELSE
454 IF(dtx(i)>dtmin1(ity).OR.off(i)==zero.
455 . or.offg(i)<zero) GO TO 75
456 ENDIF
457 off(i) = zero
458 icount=icount+1
459 list(icount)=i
460 75 CONTINUE
461
462 DO j=1,icount
463 i = list(j)
464#include "lockon.inc"
465 WRITE(iout,*)
466 . ' -- DELETE SOLID ELEMENTS',ngl(i)
467 WRITE(istdo,*)
468 . ' -- DELETE SOLID ELEMENTS',ngl(i)
469#include "lockoff.inc"
470 idel7nok = 1
471 ENDDO
472 ELSEIF(idtmin(ity)==3.AND.(ismstr==2.OR.ismstr==12))THEN
473 icount = 0
474 DO 76 i=1,nel
475 IF(dtmin(i)/=zero) THEN
476 IF(dtx(i)>dtmin(i).OR.
477 . off(i)<one.OR.offg(i)<=zero.OR.offg(i)==two) GO TO 76
478 ELSE
479 IF(dtx(i)>dtmin1(ity).OR.
480 . off(i)<one.OR.offg(i)<=zero.OR.offg(i)==two) GO TO 76
481 ENDIF
482 offg(i) = two
483 icount=icount+1
484 list(icount)=i
485 76 CONTINUE
486
487 DO j=1,icount
488 i=list(j)
489#include "lockon.inc"
490 WRITE(iout,*)
491 . '-- CONSTANT TIME STEP FOR SOLID ELEMENT NUMBER ',ngl(i)
492 WRITE(istdo,*)
493 . '-- CONSTANT TIME STEP FOR SOLID ELEMENT NUMBER ',ngl(i)
494#include "lockoff.inc"
495 ENDDO
496 ELSEIF(idtmin(ity)==5)THEN
497 error=0
498 DO 570 i=1,nel
499 IF(dtx(i)>dtmin1(ity).OR.off(i)==zero.
500 . or.offg(i)<zero) GO TO 570
501 error=1
502 570 CONTINUE
503 IF (error==1) THEN
504 mstop = 2
505 DO i=1,nel
506 IF(dtx(i)>dtmin1(ity).OR.off(i)==zero.
507 . or.offg(i)<zero) cycle
508#include "lockon.inc"
509 WRITE(iout,*)
510 . ' **ERROR : TIME STEP LESS OR EQUAL DTMIN FOR SOLID ELEMENT',
511 . ngl(i)
512#include "lockoff.inc"
513 END DO
514#include "lockon.inc"
515 WRITE(istdo,*)
516 . ' **ERROR : TIME STEP LESS OR EQUAL DTMIN FOR SOLID ELEMENTS'
517#include "lockoff.inc"
518 ENDIF
519 ENDIF
520 END IF ! IF(IMCONV==1)
521C
522 IF(idtmins == 2 .AND. jsms /= 0)THEN
523 IF(ismstr==1.OR.
524 + ((ismstr==2.OR.ismstr==12).AND.idtmin(1)==3))THEN
525 mt = mat(1)
526 DO i=1,nel
527 rho0(i)=pm(1,mt)
528 END DO
529 DO i=1,nel
530 sti(i) = zero
531 IF(off(i)==zero.OR.offg(i)<zero) cycle
532 IF(ipm(251,mt)/=0) cycle
533C 1st compute stiffness
534 IF(n2d==0) THEN
535 tidt = one/dty(i)
536 IF(offg(i)>one)THEN
537 trho = rho0(i) * tidt
538 tvol = vol0(i) * tidt
539 ELSE
540 trho = rho(i) * tidt
541 tvol = vol(i) * tidt
542 END IF
543 sti(i) = trho * tvol
544C STI will be changed to 2*STI/NNE in SxCUMU
545C [ STI(I) = FOURTH * TRHO * TVOL for 8node bricks ]
546 ELSE
547 tidt = one/dty(i)
548 trho = rho(i) * tidt
549 taire = aire(i) * tidt
550 sti(i) = half * trho * taire
551 ENDIF
552c
553c dmels = 2*dmels !!
554c w^2 < 2k / (m+dmels+dmels/...) < 2k / (m+dmels)
555c dt = 2/w = sqrt( 2*(m+dmels)/k), k=2m/dty^2
556 dmels(i)=max(dmels(i),
557 . two*mssa(i)*((dtmins/(dtfacs*dty(i)))**2 - one))
558 dtx(i)=dtfacs*sqrt(one+dmels(i)/(two*mssa(i)))*dty(i)
559
560 END DO
561 ELSE
562 DO i=1,nel
563 sti(i) = zero
564 IF(off(i)==zero.OR.offg(i)<zero) cycle
565 IF(ipm(251,mt)/=0) cycle
566C 1st compute stiffness
567 IF(n2d==0) THEN
568 tidt = one/dty(i)
569 trho = rho(i) * tidt
570 tvol = vol(i) * tidt
571 sti(i) = trho * tvol
572C STI will be changed to 2*STI/NNE in SxCUMU
573C [ STI(I) = FOURTH * TRHO * TVOL for 8node bricks ]
574 ELSE
575 tidt = one/dty(i)
576 trho = rho(i) * tidt
577 taire = aire(i) * tidt
578 sti(i) = half * trho * taire
579 ENDIF
580c
581c dmels = 2*dmels !!
582c w^2 < 2k / (m+dmels+dmels/...) < 2k / (m+dmels)
583c dt = 2/w = sqrt( 2*(m+dmels)/k), k=2m/dty^2
584 dmels(i)=max(dmels(i),
585 . two*mssa(i)*((dtmins/(dtfacs*dty(i)))**2 - one))
586 dtx(i)=dtfacs*sqrt(one+dmels(i)/(two*mssa(i)))*dty(i)
587 END DO
588 END IF
589C
590 IF(jthe > 0)THEN
591 mt = mat(1)
592 DO i=1,nel
593 IF(off(i)==zero.OR.offg(i)<zero) cycle
594 sph = pm(69,mt)
595 ak1 = pm(75,mt)
596 bk1 = pm(76,mt)
597 ak2 = pm(77,mt)
598 bk2 = pm(78,mt)
599 tli = pm(80,mt)
600 IF(temp(i)<tli)THEN
601 akk=ak1+bk1*temp(i)
602 ELSE
603 akk=ak2+bk2*temp(i)
604 ENDIF
605 IF(jtur/=0)THEN
606 xmu = rho(i)*pm(24,mt)
607 tmu = pm(81,mt)
608 rpr = pm(95,mt)
609 atu=rpr*tmu*rk(i)*rk(i)/(max(em15,re(i)*vol(i))*xmu)
610 akk=akk*(1.+atu)
611 ENDIF
612 dtx(i) = min(dtx(i),dtfacs*half*deltax(i)*deltax(i)*
613 . sph/max(akk,em20))
614 ENDDO
615 END IF
616 END IF
617C
618 IF(nodadt==0 .OR. (idtmins == 2 .AND. jsms /= 0) .OR. ale_or_euler==1)THEN
619 id_min=0
620 DO i=1,nel
621 IF(dtx(i)>dt2t .OR. off(i)<=zero .OR. offg(i)<=zero) cycle
622 ! nelts et itypts remplaces par neltst et itypst
623 IF(vol(i)<=zero)cycle
624 dt2t = dtx(i)
625 neltst = ngl(i)
626 ityptst = ity
627 id_min = i
628 ENDDO
629 IF(dt2t<dtmin1(102).AND.id_min>0)THEN
630 tstop = tt
631 print *, "DT=",dt2t,dtmin1(102)
632#include "lockon.inc"
633 WRITE(iout,*) ' **ERROR : TIME STEP LESS OR EQUAL DTMIN FOR ALE/EULER CELL',ngl(id_min)
634 WRITE(istdo,*)' **ERROR : TIME STEP LESS OR EQUAL DTMIN FOR ALE/EULER CELL',ngl(id_min)
635#include "lockoff.inc"
636 ENDIF
637 ENDIF
638C
639C--------------------------
640C THERMAL TIME STEP
641C--------------------------
642 IF (jthe < 0 .AND. glob_therm%IDT_THERM == 1)THEN
643 mt = mat(1)
644 sph = pm(69,mt)
645 ak1 = pm(75,mt)
646 bk1 = pm(76,mt)
647 ak2 = pm(77,mt)
648 bk2 = pm(78,mt)
649 tli = pm(80,mt)
650 DO i=1,nel
651 IF(temp(i)<tli)THEN
652 akk=ak1+bk1*temp(i)
653 ELSE
654 akk=ak2+bk2*temp(i)
655 ENDIF
656 IF(jtur/=0)THEN
657 xmu = rho(i)*pm(24,mt)
658 tmu = pm(81,mt)
659 rpr = pm(95,mt)
660 atu = rpr*tmu*rk(i)*rk(i)/(max(em15,re(i)*vol(i))*xmu)
661 akk = akk*(1.+atu)
662 ENDIF
663 akk = akk * glob_therm%THEACCFACT
664 dt = glob_therm%DTFACTHERM*half*deltax(i)*deltax(i)*sph/max(akk,em20)
665 IF(dt< glob_therm%DT_THERM.AND.off(i)>zero.AND.offg(i)>zero) glob_therm%DT_THERM = dt
666 conde(i) = four*vol(i)*akk/(deltax(i)*deltax(i))
667 conde(i) = conde(i)*off(i)
668 ENDDO
669 ENDIF
670C
671C--------------------------
672C ELEMENT TIME STEP
673C--------------------------
674 !here scale factor already applied
675 IF (g_dt>0)THEN
676 DO i=1,nel
677 dtel(i) = dtx(i)
678 ENDDO
679 ENDIF
680C--------------------------
681
682 RETURN
683 END
dtel
subroutine dtel(ssp, pm, geo, pid, mat, rho0, vis, deltax, aire, vol, dtx)
Definition dtel.F:46
m46law
subroutine m46law(lft, llt, nft, mtn, pm, off, sig, eint, rho, qold, vol, uvar, bufmat, stifn, mat, ngl, nuvar, dt2t, neltst, rho0, ityptst, offg, jlag, geo, pid, ssp, aire, voln, vd2, deltax, vis, d1, d2, d3, pnew, psh, q, ssp_eq, wxx, wyy, wzz, ipm, mssa, dmels, dvol, sold1, sold2, sold3, sold4, sold5, sold6, conde, vol_avg, dtel, g_dt, nel, d4, d5, d6, rhoref, rhosp, ismstr, ity, jsms, jtur, jthe, npg, svis, glob_therm)
Definition m46law.F:52
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
mmain
subroutine mmain(pm, elbuf_str, ix, nix, x, geo, iparg, nel, skew, bufmat, ipart, ipartel, nummat, matparam, imat, ipm, ngl, pid, npf, tf, mfxx, mfxy, mfxz, mfyx, mfyy, mfyz, mfzx, mfzy, mfzz, rx, ry, rz, sx, sy, sz, gama, voln, dvol, s1, s2, s3, s4, s5, s6, dxx, dyy, dzz, d4, d5, d6, wxx, wyy, wzz)
Definition mmain.F:43
mqviscb
subroutine mqviscb(pm, off, rho, rk, temp, ssp, re, sti, dt2t, neltst, ityptst, aire, offg, geo, pid, vol, vd2, deltax, vis, d1, d2, d3, pnew, psh, mat, ngl, qvis, ssp_eq, vol0, mssa, dmels, igeo, facq0, conde, dtel, g_dt, ipm, rhoref, rhosp, nel, ity, ismstr, jtur, jthe, jsms, npg, glob_therm)
Definition mqviscb.F:56
ale_mod
Definition ale_mod.F:149
ale_mod::ale
type(ale_) ale
Definition ale_mod.F:249
i22bufbric_mod
Definition cut-cell-search_mod.F:95
i22tri_mod
Definition cut-cell-search_mod.F:212