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,ICOUNT,LIST(MVSIZ),ERROR, ALE_OR_EULER, ID_MIN
112 INTEGER :: JALE_FROM_MAT, JALE_FROM_PROP
113 my_real :: DD(MVSIZ), AL(MVSIZ),
114 . DTX(MVSIZ), AD(MVSIZ), QX(MVSIZ), CX(MVSIZ), RHO0(MVSIZ), NRHO(MVSIZ),
115 . DTY(MVSIZ), QA, QB, VISI, FACQ,QAA,
116 . CNS1, CNS2, SPH, AK1, BK1, AK2, BK2, TLI, AKK, XMU, TMU, RPR,
117 . atu, dt, nui
118 my_real :: tidt,tvol,trho,taire, dtmin(mvsiz),facpg
119 my_real :: cns1_0,cns2_0,qaa_0
120#ifdef MYREAL8
121 my_real, PARAMETER :: real_three = 3.0d0
122 my_real, PARAMETER :: real_one = 1.0d0
123#else
124 my_real, PARAMETER :: real_three = 3.0
125 my_real, PARAMETER :: real_one = 1.0
126#endif
127
128C=======================================================================
129
130 ! by default don't apply /DT/NODA/* to ALE/EULER cells ; unless if /DT/NODA/ALEON is enabled (hidden card introduced for backward compatibility)
131 mt = mat(1)
132 ale_or_euler = 0
133 jale_from_mat = nint(pm(72,mt)) !ALE/EULER framework defined from /ALE/MAT or /EULER/MAT
134 jale_from_prop = igeo(62,pid(1)) !ALE/EULER framework defined from /PROP/SOLID from 3D and 2d solid elems.
135 IF(jale_from_mat == 1 .OR. jale_from_mat == 2) ale_or_euler = 1
136 IF(jale_from_prop == 1 .OR. jale_from_prop == 2)ale_or_euler = 1
137 IF(ale%GLOBAL%I_DT_NODA_ALE_ON==1) ale_or_euler = 0
138
139 IF(impl==zero)THEN
140 DO i=1,nel
141 dd(i)=-d1(i)-d2(i)-d3(i)
142 ad(i)=zero
143 al(i)=zero
144 cx(i)=ssp(i) + sqrt(vd2(i))
145 ENDDO
146 IF((impl_s>0.AND.idyna==0).OR.ismdisp>0)THEN
147 visi=zero
148 facq=zero
149 ELSE
150 visi=one
151 facq=facq0
152 ENDIF
153 IF(imconv<0) THEN
154 DO i=1,nel
155 vol(i)= abs(vol(i))
156 ENDDO
157 ENDIF
158 ELSE
159 DO i=1,nel
160 dd(i)=-d1(i)-d2(i)-d3(i)
161 ad(i)=zero
162 al(i)=zero
163 cx(i)=sqrt(vd2(i))
164 ENDDO
165 visi=zero
166 facq=zero
167 ENDIF
168C
169 IF(n2d>0) THEN
170 facpg=one
171 DO i=1,nel
172 IF(off(i)==one)THEN
173 al(i)=sqrt(aire(i))
174 ad(i)= max(zero,dd(i))
175 ENDIF
176 ENDDO
177 ELSE
178 facpg = npg
179 IF (facpg>one) facpg=facpg**(real_one/real_three)
180 DO i=1,nel
181 IF(off(i)==one)THEN
182 IF(vol(i) > zero)THEN
183 al(i) = vol(i)**(real_one/real_three)
184 ELSE
185 al(i)=zero
186 END IF
187 ad(i)= max(zero,dd(i))
188 ENDIF
189 ENDDO
190 ENDIF
191C
192 !Numerical viscosity is managed inside sigeps51.f for law51
193 !qa, qb should be 0 here otherwise viscosity will be twice.
194 IF(.NOT.(idtmins==2.AND.jsms/=0))THEN
195 mt = mat(1)
196 DO i=1,nel
197 rho0(i)=pm(1,mt)
198 nrho(i)=sqrt(rhoref(i)*rho0(i))
199 ENDDO
200 qa =facq*geo(14,pid(1))
201 qb =facq*geo(15,pid(1))
202 cns1_0=facpg*geo(16,pid(1))
203 cns2_0=facpg*geo(17,pid(1))
204 psh(1:nel)=pm(88,mt)
205 pnew(1:nel)=zero
206 qaa_0 = qa*qa
207 dtmin(1:nel) = geo(172,pid(1))
208 DO i=1,nel
209 cns1=cns1_0*al(i)*nrho(i)*ssp(i)*off(i)
210 cns2=cns2_0*al(i)*nrho(i)*ssp(i)*off(i)
211C
212C warns : Stability vs artificial viscosity should use
213C RHO(I)*AL(I)*(QAA*AL(I)+QB*SSP(I)) / MAX(EM20,RHOREF(I)*DELTAX(I))
214C Same for VIS(I) since SOUNDSPEED is generally computed wrt RHO0, not RHO ../..
215 qaa = qaa_0*ad(i)
216 qx(i)=qb*ssp(i)+al(i) * qaa
217 . + visi*two*vis(i) / max(em20,rho(i)*deltax(i))
218 . + (cns1+visi*cns2) / max(em20,rhoref(i)*deltax(i))
219 qvis(i)=rho(i)*ad(i)*al(i)*(qaa*al(i)+qb*ssp(i))
220 ENDDO
221 ELSE
222 mt = mat(1)
223 DO i=1,nel
224 rho0(i)=pm(1,mt)
225 nrho(i)=sqrt(rhoref(i)*rho0(i))
226 ENDDO
227
228 qa =facq*geo(14,pid(1))
229 qb =facq*geo(15,pid(1))
230 psh(1:nel)=pm(88,mt)
231 pnew(1:nel)=zero
232 qaa_0 = qa*qa
233 dtmin(1:nel) = geo(172,pid(1))
234
235 IF(geo(16,pid(1)) >= zero)THEN
236 cns1_0=facpg*geo(16,pid(1))
237 cns2_0=facpg*geo(17,pid(1))
238 DO i=1,nel
239 cns1=cns1_0*al(i)*nrho(i)*ssp(i)*off(i)
240 cns2=cns2_0*al(i)*nrho(i)*ssp(i)*off(i)
241 qaa = qaa_0*ad(i)
242 nui =rho(i)*al(i)*(qaa*al(i)+qb*ssp(i))
243 qvis(i)=nui*ad(i)
244 qx(i) =(nui + cns1+visi*(two*vis(i)+cns2))/max(em20,rhoref(i)*deltax(i))
245 ENDDO
246 ELSE
247 cns1_0=abs(geo(16,pid(1)))
248 cns2_0=abs(geo(17,pid(1)))
249 DO i=1,nel
250 cns1=cns1_0*nrho(i)*ssp(i)**2*off(i)
251 cns2=cns2_0*nrho(i)*ssp(i)**2*off(i)
252 qaa = qaa_0*ad(i)
253 nui =rho(i)*al(i)*(qaa*al(i)+qb*ssp(i))
254 qvis(i)=nui*ad(i)
255 qx(i) =(nui + cns1+visi*(two*vis(i)+cns2))/max(em20,rhoref(i)*deltax(i))
256 ENDDO
257 END IF
258 ENDIF
259C
260 DO i=1,nel
261 ssp_eq(i) = max(em20,qx(i)+sqrt(qx(i)*qx(i)+cx(i)*cx(i)))
262 dtx(i) = deltax(i) / ssp_eq(i)
263 ENDDO
264C
265 !This would output CFL value without scale factor.
266 !IF (G_DT>0)THEN
267 ! DO I=1,NEL
268 ! DTEL(I) = DTX(I)
269 ! ENDDO
270 !ENDIF
271C
272 IF(int22 > 0)THEN
273 !check if at least one brick has time step smaller than cinematic one => elementary time step.
274 DO i=1,nel
275 IF(dtx(i)<dt22_min)THEN
276 idt_int22 = 0
277 !NELTST = NGL(I)
278 !ITYPTST = ITY
279 EXIT ! just set IDT_INT22
280 ENDIF
281 ENDDO
282 ENDIF
283C
284 IF(jthe > 0)THEN
285 mt = mat(1)
286 sph = pm(69,mt)
287 ak1 = pm(75,mt)
288 bk1 = pm(76,mt)
289 ak2 = pm(77,mt)
290 bk2 = pm(78,mt)
291 tli = pm(80,mt)
292 DO i=1,nel
293 IF(temp(i)<tli)THEN
294 akk=ak1+bk1*temp(i)
295 ELSE
296 akk=ak2+bk2*temp(i)
297 ENDIF
298 IF(jtur/=0)THEN
299 xmu = rho(i)*pm(24,mt)
300 tmu = pm(81,mt)
301 rpr = pm(95,mt)
302 atu=rpr*tmu*rk(i)*rk(i)/(max(em15,re(i)*vol(i))*xmu)
303 akk=akk*(1.+atu)
304 ENDIF
305 dtx(i) = min(dtx(i),half*deltax(i)*deltax(i)*sph/max(akk,em20))
306 ENDDO
307 ENDIF
308C
309 IF(.NOT.(idtmins==2.AND.jsms/=0))THEN
310 IF(kdtsmstr==1.AND.ismstr==1.OR.
311 . ((ismstr==2.OR.ismstr==12).AND.idtmin(1)==3))THEN
312 ! Kdtsmstr == 1 in version 5, by default.
313 mt = mat(1)
314 DO i=1,nel
315 rho0(i)=pm(1,mt)
316 END DO
317 DO 50 i=1,nel
318 sti(i) = zero
319 IF(off(i)==zero.OR.offg(i)<zero) GO TO 50
320 IF(ipm(251,mt)/=0) GO TO 50
321 IF(n2d==0) THEN
322 tidt = one/dtx(i)
323 IF(offg(i)>one)THEN
324 trho = rho0(i) * tidt
325 tvol = vol0(i) * tidt
326 ELSE
327 trho = rho(i) * tidt
328 tvol = vol(i) * tidt
329 END IF
330 ! STI will be changed to 2*STI/NNE in SxCUMU
331 ! [ STI(I) = FOURTH * TRHO * TVOL for 8node bricks ]
332 sti(i) = trho * tvol
333 ELSE
334 ! small strain is not available in 2D analysis
335 tidt = one/dtx(i)
336 trho = rho(i) * tidt
337 taire = aire(i) * tidt
338 sti(i) = half * trho * taire
339 ENDIF
340C dt2 replaced by dt2t
341 50 CONTINUE
342 IF(ale_or_euler==0)THEN
343 DO i=1,nel
344 dtx(i)= dtfac1(ity)*dtx(i)
345 ENDDO
346 ELSE
347 DO i=1,nel
348 dtx(i)= dtfac1(102)*dtx(i)
349 ENDDO
350 ENDIF
351 !/DT/NODA & /DT/NODA/CST HAS NO EFFECT WITH ALE/EULER
352 IF(ale_or_euler==1 .OR. nodadt==0)THEN
353 DO i=1,nel
354 IF(vol(i)>zero .AND. (off(i)/=zero.AND.offg(i)>=zero))dt2t= min(dtx(i),dt2t)
355 ENDDO
356 ENDIF
357
358 ELSE
359 IF(ismstr == 11) THEN
360 IF(n2d == 0) THEN
361 DO i=1,nel
362 sti(i) = zero
363 IF(off(i)==zero.OR.offg(i)<zero) cycle
364 IF(ipm(251,mt)/=0) cycle
365 tidt = one/dtx(i)
366 trho = rho0(i) * tidt
367 tvol = vol0(i) * tidt
368 sti(i) = trho * tvol
369 ENDDO
370 ELSE
371 DO i=1,nel
372 sti(i) = zero
373 IF(off(i)==zero.OR.offg(i)<zero) cycle
374 IF(ipm(251,mt)/=0) cycle
375 tidt = one/dtx(i)
376 trho = rho(i) * tidt
377 taire = aire(i) * tidt
378 sti(i) = half * trho * taire
379 ENDDO
380 ENDIF ! N2D
381 ELSE
382 DO 60 i=1,nel
383 sti(i) = zero
384 IF(off(i)==zero.OR.offg(i)<zero) GO TO 60
385 IF(ipm(251,mt)/=0) GO TO 60
386 IF(n2d==0) THEN
387 tidt = one/dtx(i)
388 trho = rho(i) * tidt
389 tvol = vol(i) * tidt
390 sti(i) = trho * tvol
391C STI will be changed to 2*STI/NNE in SxCUMU
392C [ STI(I) = FOURTH * TRHO * TVOL for 8node bricks ]
393 ELSE
394 tidt = one/dtx(i)
395 trho = rho(i) * tidt
396 taire = aire(i) * tidt
397 sti(i) = half * trho * taire
398 ENDIF
399C dt2 replaced by dt2t
400 60 CONTINUE
401 ENDIF ! ISMSTR == 11
402 IF(ale_or_euler==0)THEN
403 DO i=1,nel
404 dtx(i)= dtfac1(ity)*dtx(i)
405 ENDDO
406 ELSE
407 DO i=1,nel
408 dtx(i)= dtfac1(102)*dtx(i)
409 ENDDO
410 ENDIF
411 !/DT/NODA & /DT/NODA/CST HAS NO EFFECT WITH ALE/EULER
412 IF(ale_or_euler==1 .OR. nodadt==0)THEN
413 DO i=1,nel
414 IF(vol(i)>zero .AND. (off(i)/=zero.AND.offg(i)>=zero))dt2t= min(dtx(i),dt2t)
415 ENDDO
416 ENDIF
417 END IF
418C
419 ELSE ! IDTMINS == 2 .AND. JSMS /= 0
420 DO i=1,nel
421 dty(i)= dtx(i)
422 dtx(i)= dtfac1(ity)*dtx(i)
423 END DO
424 END IF
425C----
426 IF(imconv==1)THEN
427 IF(idtmin(ity)==1)THEN
428 error=0
429 DO 70 i=1,nel
430 IF(dtx(i)>dtmin1(ity).OR.off(i)==zero.
431 . or.offg(i)<zero) GO TO 70
432 error=1
433 70 CONTINUE
434
435 IF (error==1) THEN
436 tstop = tt
437 DO i=1,nel
438 IF(dtx(i)>dtmin1(ity).OR.off(i)==zero.
439 . or.offg(i)<zero) cycle
440#include "lockon.inc"
441 WRITE(iout,*)
442 . ' **ERROR : TIME STEP LESS OR EQUAL DTMIN FOR SOLID ELEMENT',
443 . ngl(i)
444#include "lockoff.inc"
445 END DO
446#include "lockon.inc"
447 WRITE(istdo,*)
448 . ' **ERROR : TIME STEP LESS OR EQUAL DTMIN FOR SOLID ELEMENTS'
449#include "lockoff.inc"
450 ENDIF
451 ELSEIF(idtmin(ity)==2)THEN
452 icount=0
453 DO 75 i=1,nel
454 IF(dtmin(i)/=zero) THEN
455 IF(dtx(i)>dtmin(i).OR.off(i)==zero.
456 . or.offg(i)<zero) GO TO 75
457 ELSE
458 IF(dtx(i)>dtmin1(ity).OR.off(i)==zero.
459 . or.offg(i)<zero) GO TO 75
460 ENDIF
461 off(i) = zero
462 icount=icount+1
463 list(icount)=i
464 75 CONTINUE
465
466 DO j=1,icount
467 i = list(j)
468#include "lockon.inc"
469 WRITE(iout,*)
470 . ' -- DELETE SOLID ELEMENTS',ngl(i)
471 WRITE(istdo,*)
472 . ' -- DELETE SOLID ELEMENTS',ngl(i)
473#include "lockoff.inc"
474 idel7nok = 1
475 ENDDO
476 ELSEIF(idtmin(ity)==3.AND.(ismstr==2.OR.ismstr==12))THEN
477 icount = 0
478 DO 76 i=1,nel
479 IF(dtmin(i)/=zero) THEN
480 IF(dtx(i)>dtmin(i).OR.
481 . off(i)<one.OR.offg(i)<=zero.OR.offg(i)==two) GO TO 76
482 ELSE
483 IF(dtx(i)>dtmin1(ity).OR.
484 . off(i)<one.OR.offg(i)<=zero.OR.offg(i)==two) GO TO 76
485 ENDIF
486 offg(i) = two
487 icount=icount+1
488 list(icount)=i
489 76 CONTINUE
490
491 DO j=1,icount
492 i=list(j)
493#include "lockon.inc"
494 WRITE(iout,*)
495 . '-- CONSTANT TIME STEP FOR SOLID ELEMENT NUMBER ',ngl(i)
496 WRITE(istdo,*)
497 . '-- CONSTANT TIME STEP FOR SOLID ELEMENT NUMBER ',ngl(i)
498#include "lockoff.inc"
499 ENDDO
500 ELSEIF(idtmin(ity)==5)THEN
501 error=0
502 DO 570 i=1,nel
503 IF(dtx(i)>dtmin1(ity).OR.off(i)==zero.
504 . or.offg(i)<zero) GO TO 570
505 error=1
506 570 CONTINUE
507 IF (error==1) THEN
508 mstop = 2
509 DO i=1,nel
510 IF(dtx(i)>dtmin1(ity).OR.off(i)==zero.
511 . or.offg(i)<zero) cycle
512#include "lockon.inc"
513 WRITE(iout,*)
514 . ' **ERROR : TIME STEP LESS OR EQUAL DTMIN FOR SOLID ELEMENT',
515 . ngl(i)
516#include "lockoff.inc"
517 END DO
518#include "lockon.inc"
519 WRITE(istdo,*)
520 . ' **ERROR : TIME STEP LESS OR EQUAL DTMIN FOR SOLID ELEMENTS'
521#include "lockoff.inc"
522 ENDIF
523 ENDIF
524 END IF ! IF(IMCONV==1)
525C
526 IF(idtmins == 2 .AND. jsms /= 0)THEN
527 IF(ismstr==1.OR.
528 + ((ismstr==2.OR.ismstr==12).AND.idtmin(1)==3))THEN
529 mt = mat(1)
530 DO i=1,nel
531 rho0(i)=pm(1,mt)
532 END DO
533 DO i=1,nel
534 sti(i) = zero
535 IF(off(i)==zero.OR.offg(i)<zero) cycle
536 IF(ipm(251,mt)/=0) cycle
537C 1st compute stiffness
538 IF(n2d==0) THEN
539 tidt = one/dty(i)
540 IF(offg(i)>one)THEN
541 trho = rho0(i) * tidt
542 tvol = vol0(i) * tidt
543 ELSE
544 trho = rho(i) * tidt
545 tvol = vol(i) * tidt
546 END IF
547 sti(i) = trho * tvol
548C STI will be changed to 2*STI/NNE in SxCUMU
549C [ STI(I) = FOURTH * TRHO * TVOL for 8node bricks ]
550 ELSE
551 tidt = one/dty(i)
552 trho = rho(i) * tidt
553 taire = aire(i) * tidt
554 sti(i) = half * trho * taire
555 ENDIF
556c
557c dmels = 2*dmels !!
558c w^2 < 2k / (m+dmels+dmels/...) < 2k / (m+dmels)
559c dt = 2/w = sqrt( 2*(m+dmels)/k), k=2m/dty^2
560 dmels(i)=max(dmels(i),
561 . two*mssa(i)*((dtmins/(dtfacs*dty(i)))**2 - one))
562 dtx(i)=dtfacs*sqrt(one+dmels(i)/(two*mssa(i)))*dty(i)
563
564 END DO
565 ELSE
566 DO i=1,nel
567 sti(i) = zero
568 IF(off(i)==zero.OR.offg(i)<zero) cycle
569 IF(ipm(251,mt)/=0) cycle
570C 1st compute stiffness
571 IF(n2d==0) THEN
572 tidt = one/dty(i)
573 trho = rho(i) * tidt
574 tvol = vol(i) * tidt
575 sti(i) = trho * tvol
576C STI will be changed to 2*STI/NNE in SxCUMU
577C [ STI(I) = FOURTH * TRHO * TVOL for 8node bricks ]
578 ELSE
579 tidt = one/dty(i)
580 trho = rho(i) * tidt
581 taire = aire(i) * tidt
582 sti(i) = half * trho * taire
583 ENDIF
584c
585c dmels = 2*dmels !!
586c w^2 < 2k / (m+dmels+dmels/...) < 2k / (m+dmels)
587c dt = 2/w = sqrt( 2*(m+dmels)/k), k=2m/dty^2
588 dmels(i)=max(dmels(i),
589 . two*mssa(i)*((dtmins/(dtfacs*dty(i)))**2 - one))
590 dtx(i)=dtfacs*sqrt(one+dmels(i)/(two*mssa(i)))*dty(i)
591 END DO
592 END IF
593C
594 IF(jthe > 0)THEN
595 mt = mat(1)
596 DO i=1,nel
597 IF(off(i)==zero.OR.offg(i)<zero) cycle
598 sph = pm(69,mt)
599 ak1 = pm(75,mt)
600 bk1 = pm(76,mt)
601 ak2 = pm(77,mt)
602 bk2 = pm(78,mt)
603 tli = pm(80,mt)
604 IF(temp(i)<tli)THEN
605 akk=ak1+bk1*temp(i)
606 ELSE
607 akk=ak2+bk2*temp(i)
608 ENDIF
609 IF(jtur/=0)THEN
610 xmu = rho(i)*pm(24,mt)
611 tmu = pm(81,mt)
612 rpr = pm(95,mt)
613 atu=rpr*tmu*rk(i)*rk(i)/(max(em15,re(i)*vol(i))*xmu)
614 akk=akk*(1.+atu)
615 ENDIF
616 dtx(i) = min(dtx(i),dtfacs*half*deltax(i)*deltax(i)*
617 . sph/max(akk,em20))
618 ENDDO
619 END IF
620 END IF
621C
622 IF(nodadt==0 .OR. (idtmins == 2 .AND. jsms /= 0) .OR. ale_or_euler==1)THEN
623 id_min=0
624 DO i=1,nel
625 IF(dtx(i)>dt2t .OR. off(i)<=zero .OR. offg(i)<=zero) cycle
626 ! nelts and itypts replaced by neltst and itypst
627 IF(vol(i)<=zero)cycle
628 dt2t = dtx(i)
629 neltst = ngl(i)
630 ityptst = ity
631 id_min = i
632 ENDDO
633 IF(dt2t<dtmin1(102).AND.id_min>0)THEN
634 tstop = tt
635 print *, "DT=",dt2t,dtmin1(102)
636#include "lockon.inc"
637 WRITE(iout,*) ' **ERROR : TIME STEP LESS OR EQUAL DTMIN FOR ALE/EULER CELL',ngl(id_min)
638 WRITE(istdo,*)' **ERROR : TIME STEP LESS OR EQUAL DTMIN FOR ALE/EULER CELL',ngl(id_min)
639#include "lockoff.inc"
640 ENDIF
641 ENDIF
642C
643C--------------------------
644C THERMAL TIME STEP
645C--------------------------
646 IF (jthe < 0 .AND. glob_therm%IDT_THERM == 1)THEN
647 mt = mat(1)
648 sph = pm(69,mt)
649 ak1 = pm(75,mt)
650 bk1 = pm(76,mt)
651 ak2 = pm(77,mt)
652 bk2 = pm(78,mt)
653 tli = pm(80,mt)
654 DO i=1,nel
655 IF(temp(i)<tli)THEN
656 akk=ak1+bk1*temp(i)
657 ELSE
658 akk=ak2+bk2*temp(i)
659 ENDIF
660 IF(jtur/=0)THEN
661 xmu = rho(i)*pm(24,mt)
662 tmu = pm(81,mt)
663 rpr = pm(95,mt)
664 atu = rpr*tmu*rk(i)*rk(i)/(max(em15,re(i)*vol(i))*xmu)
665 akk = akk*(1.+atu)
666 ENDIF
667 akk = akk * glob_therm%THEACCFACT
668 dt = glob_therm%DTFACTHERM*half*deltax(i)*deltax(i)*sph/max(akk,em20)
669 IF(dt< glob_therm%DT_THERM.AND.off(i)>zero.AND.offg(i)>zero) glob_therm%DT_THERM = dt
670 conde(i) = four*vol(i)*akk/(deltax(i)*deltax(i))
671 conde(i) = conde(i)*off(i)
672 ENDDO
673 ENDIF
674C
675C--------------------------
676C ELEMENT TIME STEP
677C--------------------------
678 !here scale factor already applied
679 IF (g_dt>0)THEN
680 DO i=1,nel
681 dtel(i) = dtx(i)
682 ENDDO
683 ENDIF
684C--------------------------
685
686 RETURN
687 END
dtel
subroutine dtel(ssp, pm, geo, pid, mat, rho0, vis, deltax, aire, vol, dtx)
Definition dtel.F:46
m22law
subroutine m22law(pm, off, sig, eint, rho, qold, epxe, epd, vol, pdam, stifn, dt2t, neltst, ityptst, offg, geo, pid, amu, vol_avg, mumax, mat, ngl, ssp, dvol, aire, vnew, vd2, deltax, vis, d1, d2, d3, d4, d5, d6, pnew, psh, qnew, ssp_eq, sold1, sold2, sold3, sold4, sold5, sold6, ipla, sigy, defp, dpla, mssa, dmels, conde, dtel, g_dt, nel, ipm, rhoref, rhosp, nft, jsph, ity, jtur, jthe, ismstr, jsms, npg, glob_therm, numgeo, igeo)
Definition m22law.F:53
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
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:150
ale_mod::ale
type(ale_) ale
Definition ale_mod.F:253
i22bufbric_mod
Definition cut-cell-search_mod.F:95
i22tri_mod
Definition cut-cell-search_mod.F:212