OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mdtsph.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!|| mdtsph ../engine/source/materials/mat_share/mdtsph.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!|| mmain ../engine/source/materials/mat_share/mmain.F90
33!|| mulaw ../engine/source/materials/mat_share/mulaw.F90
34!|| sboltlaw ../engine/source/elements/solid/solide/sboltlaw.F
35!|| usermat_solid ../engine/source/materials/mat_share/usermat_solid.F
36!||====================================================================
37 SUBROUTINE mdtsph(
38 1 PM, OFF, RHO, RK,
39 2 T, RE, STI, DT2T,
40 3 NELTST, ITYPTST, OFFG, GEO,
41 4 PID, MUMAX, SSP, VOL,
42 5 VD2, DELTAX, VIS, D1,
43 6 D2, D3, PNEW, PSH,
44 7 MAT, NGL, QVIS, SSP_EQ,
45 8 G_DT, DTSPH, NEL, ITY,
46 9 JTUR, JTHE)
47C-----------------------------------------------
48C I m p l i c i t T y p e s
49C-----------------------------------------------
50#include "implicit_f.inc"
51#include "comlock.inc"
52C-----------------------------------------------
53C G l o b a l P a r a m e t e r s
54C-----------------------------------------------
55#include "mvsiz_p.inc"
56C-----------------------------------------------
57C C o m m o n B l o c k s
58C-----------------------------------------------
59#include "com08_c.inc"
60#include "scr02_c.inc"
61#include "scr18_c.inc"
62#include "param_c.inc"
63#include "cong1_c.inc"
64#include "units_c.inc"
65#include "scr07_c.inc"
66C-----------------------------------------------
67C D u m m y A r g u m e n t s
68C-----------------------------------------------
69 INTEGER, INTENT(IN) :: NEL
70 INTEGER, INTENT(IN) :: ITY
71 INTEGER, INTENT(IN) :: JTUR
72 INTEGER, INTENT(IN) :: JTHE
73 INTEGER NELTST,ITYPTST,PID(*),MAT(*), NGL(*)
74 my_real :: DT2T
75
76 my_real
77 . PM(NPROPM,*), OFF(*), RHO(*), RK(*), T(*),
78 . re(*),sti(*),offg(*),geo(npropg,*),mumax(*),
79 . vol(*), vd2(*), deltax(*), ssp(*), vis(*),
80 . psh(*), pnew(*),qvis(*) ,ssp_eq(*), d1(*),
81 . d2(*), d3(*)
82 my_real, INTENT(INOUT) :: dtsph(1:nel)
83 INTEGER,INTENT(IN) :: G_DT
84C-----------------------------------------------
85C L o c a l V a r i a b l e s
86C-----------------------------------------------
87 INTEGER I,J, MT
88 my_real
89 . AL(MVSIZ),DTX(MVSIZ), QX(MVSIZ), CX(MVSIZ), QXMATER(MVSIZ),
90 . QA, QB, VISI, FACQ,
91 . cns1, cns2, sph, ak1, bk1, ak2, bk2, tli, akk, xmu, tmu, rpr,
92 . atu
93C-----------------------------------------------
94C S o u r c e L i n e s
95C-----------------------------------------------
96 IF(impl==zero)THEN
97 DO i=1,nel
98 cx(i)=ssp(i)+sqrt(vd2(i))
99 ENDDO
100 visi=one
101 facq=one
102 ELSE
103 DO i=1,nel
104 cx(i)=sqrt(vd2(i))
105 ENDDO
106 visi=zero
107 facq=zero
108 ENDIF
109
110 !not a bug : only law 24 uses CNS1 & CNS2
111 !(they are not yet available with SPH).
112 DO i=1,nel
113 al(i)=zero
114 IF(off(i)<1.) cycle
115 al(i)=vol(i)**third
116 ENDDO
117
118 mt = mat(1)
119 DO i=1,nel
120 qa =facq*geo(14,pid(i))
121 qb =facq*geo(15,pid(i))
122 cns1=geo(16,pid(i))
123 cns2=geo(17,pid(i))*ssp(i)*al(i)*rho(i)
124 psh(i)=pm(88,mt)
125 pnew(i)=zero
126 qxmater(i)=cns1*ssp(i) + visi*(two*vis(i)+cns2) / max(em20,rho(i)*deltax(i))
127 qx(i)=qb*ssp(i) + qa*mumax(i) + qxmater(i)
128 qvis(i)=zero
129 ENDDO
130
131 DO i=1,nel
132 dtx(i)=deltax(i)/max(em20,qx(i)+sqrt(qx(i)*qx(i)+cx(i)*cx(i)))
133 !preparing material sound speed for nodal time step computation:
134 ssp_eq(i) = max(em20,qxmater(i)+sqrt(qxmater(i)*qxmater(i)+cx(i)*cx(i)))
135 ENDDO
136
137 IF(jthe/=0)THEN
138 mt = mat(1)
139 sph = pm(69,mt)
140 ak1 = pm(75,mt)
141 bk1 = pm(76,mt)
142 ak2 = pm(77,mt)
143 bk2 = pm(78,mt)
144 tli = pm(80,mt)
145 DO i=1,nel
146 IF(t(i)<tli)THEN
147 akk=ak1+bk1*t(i)
148 ELSE
149 akk=ak2+bk2*t(i)
150 ENDIF
151 IF(jtur/=0)THEN
152 xmu = rho(i)*pm(24,mt)
153 tmu = pm(81,mt)
154 rpr = pm(95,mt)
155 atu=rpr*tmu*rk(i)*rk(i)/(max(em15,re(i)*vol(i))*xmu)
156 akk=akk*(one + atu)
157 ENDIF
158 dtx(i) = min(dtx(i),half*deltax(i)*deltax(i)*sph/akk)
159 ENDDO
160 ENDIF
161
162 DO i=1,nel
163 sti(i) = zero
164 IF(off(i)==zero) cycle
165 sti(i) = two*rho(i) * vol(i) / (dtx(i)*dtx(i))
166 dtx(i)= dtfac1(ity)*dtx(i)
167 !dt2 remplace par dt2t
168 IF(nodadt==0)dt2t= min(dtx(i),dt2t)
169 ENDDO
170
171 IF(g_dt/=zero)THEN
172 DO i=1,nel
173 dtsph(i) = dtx(i)
174 ENDDO
175 ENDIF
176
177
178 IF(nodadt==0)THEN
179 IF(idtmin(ity)==1)THEN
180 DO 170 i=1,nel
181 IF(dtx(i)>dtmin1(ity).OR.off(i)==zero) GO TO 170
182 tstop = tt
183#include "lockon.inc"
184 WRITE(iout,*) ' **ERROR : TIME STEP LESS OR EQUAL DTMIN FOR SPH PARTICLE'
185 WRITE(istdo,*)' **ERROR : TIME STEP LESS OR EQUAL DTMIN FOR SPH PARTICLE'
186#include "lockoff.inc"
187 170 CONTINUE
188 ELSEIF(idtmin(ity)==2)THEN
189 DO 270 i=1,nel
190 IF(dtx(i)>dtmin1(ity).OR.off(i)==zero) GO TO 270
191 off(i) = zero
192#include "lockon.inc"
193 WRITE(iout,*) ' -- DELETE SPH PARTICLE',ngl(i)
194 WRITE(istdo,*)' -- DELETE SPH PARTICLE',ngl(i)
195#include "lockoff.inc"
196 270 CONTINUE
197 ELSEIF(idtmin(ity)==5)THEN
198 DO 570 i=1,nel
199 IF(dtx(i)>dtmin1(ity).OR.off(i)==0.) GO TO 570
200 mstop = 2
201#include "lockon.inc"
202 WRITE(iout,*)
203 . ' **ERROR : TIME STEP LESS OR EQUAL DTMIN FOR SPH PARTICLE'
204 WRITE(istdo,*)
205 . ' **ERROR : TIME STEP LESS OR EQUAL DTMIN FOR SPH PARTICLE'
206#include "lockoff.inc"
207 570 CONTINUE
208 ENDIF
209
210 DO i=1,nel
211 IF(dtx(i)>dt2t.OR.off(i)==zero) cycle
212 !nelts et itypts remplaces par neltst et itypst
213 neltst =ngl(i)
214 ityptst=ity
215 ENDDO
216
217 ENDIF
218
219 RETURN
220 END
dtsph
subroutine dtsph(ssp, pm, geo, pid, mat, rho0, vis, deltax, vol, dtx)
Definition dtsph.F:44
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
mdtsph
subroutine mdtsph(pm, off, rho, rk, t, re, sti, dt2t, neltst, ityptst, offg, geo, pid, mumax, ssp, vol, vd2, deltax, vis, d1, d2, d3, pnew, psh, mat, ngl, qvis, ssp_eq, g_dt, dtsph, nel, ity, jtur, jthe)
Definition mdtsph.F:47