OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mqviscb.F File Reference
#include "implicit_f.inc"
#include "comlock.inc"
#include "mvsiz_p.inc"
#include "com01_c.inc"
#include "com04_c.inc"
#include "com08_c.inc"
#include "cong1_c.inc"
#include "impl1_c.inc"
#include "param_c.inc"
#include "scr02_c.inc"
#include "scr07_c.inc"
#include "scr17_c.inc"
#include "scr18_c.inc"
#include "sms_c.inc"
#include "units_c.inc"
#include "inter22.inc"

Go to the source code of this file.

Functions/Subroutines

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)

Function/Subroutine Documentation

◆ mqviscb()

subroutine mqviscb ( pm,
off,
rho,
rk,
temp,
ssp,
re,
sti,
dt2t,
integer neltst,
integer ityptst,
aire,
offg,
geo,
integer, dimension(*) pid,
vol,
vd2,
deltax,
vis,
d1,
d2,
d3,
pnew,
psh,
integer, dimension(*) mat,
integer, dimension(*) ngl,
qvis,
ssp_eq,
vol0,
mssa,
dmels,
integer, dimension(npropgi,numgeo) igeo,
facq0,
conde,
dimension(mvsiz), intent(inout) dtel,
integer g_dt,
integer, dimension(npropmi,nummat) ipm,
dimension(mvsiz), intent(in) rhoref,
dimension(mvsiz), intent(in) rhosp,
integer, intent(in) nel,
integer, intent(in) ity,
integer, intent(in) ismstr,
integer, intent(in) jtur,
integer, intent(in) jthe,
integer, intent(in) jsms,
integer, intent(in) npg,
type (glob_therm_), intent(inout) glob_therm )

Definition at line 43 of file mqviscb.F.

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
my_real
#define my_real
Definition cppsort.cpp:32
dtel
subroutine dtel(ssp, pm, geo, pid, mat, rho0, vis, deltax, aire, vol, dtx)
Definition dtel.F:46
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
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