OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sboltlaw.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!|| sboltlaw ../engine/source/elements/solid/solide/sboltlaw.F
25!||--- called by ------------------------------------------------------
26!|| mmain ../engine/source/materials/mat_share/mmain.F90
27!||--- calls -----------------------------------------------------
28!|| mdtsph ../engine/source/materials/mat_share/mdtsph.F
29!|| mqviscb ../engine/source/materials/mat_share/mqviscb.F
30!||--- uses -----------------------------------------------------
31!|| glob_therm_mod ../common_source/modules/mat_elem/glob_therm_mod.F90
32!|| prop_param_mod ../common_source/modules/mat_elem/prop_param_mod.F90
33!||====================================================================
34 SUBROUTINE sboltlaw(
35 1 PM, SIG, MAT, D1,
36 2 D2, D3, D4, D5,
37 3 D6, NEL, RHO, BPRELD,
38 4 EINT, QOLD, VOL0, STIFN,
39 5 DT2T, NELTST, ITYPTST, OFFG,
40 6 GEO, PID, MUMAX, NGL,
41 7 SSP, DVOL, AIRE, VNEW,
42 8 VD2, DELTAX, VIS, PNEW,
43 9 PSH, QNEW, SSP_EQ, SOLD1,
44 A SOLD2, SOLD3, SOLD4, SOLD5,
45 B SOLD6, MSSA, DMELS, CONDE,
46 C AMU, VOL_AVG, DTEL, G_DT,
47 D OFF, IPM, RHOREF, RHOSP,
48 E VOL0DP, ISMSTR, JSPH, JTUR,
49 F ITY, JTHE, JSMS, NPG ,
50 G glob_therm,NUMGEO,IGEO)
51C-----------------------------------------------
52C M o d u l e s
53C-----------------------------------------------
54 use glob_therm_mod
55 USE prop_param_mod , only : n_var_igeo
56C-----------------------------------------------
57C I m p l i c i t T y p e s
58C-----------------------------------------------
59#include "implicit_f.inc"
60C-----------------------------------------------
61C G l o b a l P a r a m e t e r s
62C-----------------------------------------------
63#include "mvsiz_p.inc"
64C-----------------------------------------------
65C C o m m o n B l o c k s
66C-----------------------------------------------
67#include "com08_c.inc"
68#include "param_c.inc"
69#include "scr05_c.inc"
70C-----------------------------------------------
71C D u m m y A r g u m e n t s
72C-----------------------------------------------
73 INTEGER, INTENT(IN) :: JSMS
74 INTEGER, INTENT(IN) :: JTUR
75 INTEGER, INTENT(IN) :: ITY
76 INTEGER, INTENT(IN) :: JTHE
77 INTEGER, INTENT(IN) :: ISMSTR
78 INTEGER, INTENT(IN) :: JSPH,NPG
79 INTEGER,INTENT(IN) :: NUMGEO
80C
81 INTEGER NEL
82 my_real
83 . PM(NPROPM,*), SIG(NEL,6), RHO(MVSIZ),
84 . D1(*), D2(*), D3(*), D4(*),D5(*), D6(*),
85 . BPRELD(NEL,*), OFF(*)
86 INTEGER NELTST,ITYPTST,PID(*),G_DT
87 INTEGER MAT(*),NGL(*),IPM(NPROPMI,*)
88 my_real
89 . dt2t
90
91 my_real
92 . eint(*), qold(*),
93 . vol0(*),stifn(*), offg(*),geo(npropg,*),mumax(*)
94 my_real
95 . vnew(*), vd2(*), deltax(*), ssp(*), aire(*), vis(*),
96 . psh(*), pnew(*),qnew(*) ,ssp_eq(*), dvol(*),
97 . sold1(*), sold2(*), sold3(*), sold4(*), sold5(*), sold6(*),
98 . mssa(*), dmels(*),conde(*),amu(*),vol_avg(*),dtel(*),
99 . rhoref(*) ,rhosp(*)
100 DOUBLE PRECISION VOL0DP(*)
101 type (glob_therm_) ,intent(inout) :: glob_therm
102 integer,dimension(n_var_igeo,numgeo),intent(in) :: igeo
103C-----------------------------------------------
104C L o c a l V a r i a b l e s
105C-----------------------------------------------
106 INTEGER I, MX !, J
107 my_real
108 . reduc ,reduc1 , t1, t2 , ee, nu, ggt, treps, p, tp
109 my_real
110 . c1(mvsiz), lmbd(mvsiz), gg(mvsiz),
111 . rho0(mvsiz), rho0_1,c1_1
112C-----------------------------------------------
113 my_real
114 . facq0,
115 . e1, e2, e3, e4, e5, e6, einc,
116 . bid1, bid2, bid3, dta
117C-----------------------------------------------
118 facq0 = one
119 reduc1 = em04 !EM03 !ONE
120
121 DO i=1,nel
122 t1 = bpreld(i,1)+zep4*(bpreld(i,2)-bpreld(i,1))
123 t2 = bpreld(i,1)+zep7*(bpreld(i,2)-bpreld(i,1))
124 reduc = reduc1
125 tp = tt-t1
126 IF(tp > zero ) THEN
127 reduc = min(reduc1*(one-tp/(t2-t1))+tp/(t2-t1), one)
128 ENDIF
129C
130 mx = mat(i)
131 rho0_1 = pm(1, mx)
132 nu = pm(21, mx)
133 ee = pm(20, mx)*reduc
134 gg(i) = pm(22, mx)*reduc
135 c1_1 = pm(32, mx)
136 lmbd(i) = dt1*ee*nu/((one+nu)*(one-two*nu))
137 c1(i) = c1_1*reduc
138 rho0(i) =rho0_1
139C
140 !SSP(I)=SQRT((ONEP333*GG(I)+C1(I))/RHO0(I))
141 ssp(i)=sqrt((onep333*pm(22, mx)+c1_1)/rho0(i))
142 ENDDO
143C
144 IF (jsph==0)THEN
145 CALL mqviscb(
146 1 pm, off, rho, bid1,
147 2 bid2, ssp, bid3, stifn,
148 3 dt2t, neltst, ityptst, aire,
149 4 offg, geo, pid, vnew,
150 5 vd2, deltax, vis, d1,
151 6 d2, d3, pnew, psh,
152 7 mat, ngl, qnew, ssp_eq,
153 8 vol0, mssa, dmels, igeo,
154 9 facq0, conde, dtel, g_dt,
155 a ipm, rhoref, rhosp, nel,
156 b ity, ismstr, jtur, jthe,
157 c jsms, npg , glob_therm)
158 ELSE
159 CALL mdtsph(
160 1 pm, off, rho, bid1,
161 2 bid2, bid3, stifn, dt2t,
162 3 neltst, ityptst, offg, geo,
163 4 pid, mumax, ssp, vnew,
164 5 vd2, deltax, vis, d1,
165 6 d2, d3, pnew, psh,
166 7 mat, ngl, qnew, ssp_eq,
167 8 g_dt, dtel, nel, ity,
168 9 jtur, jthe)
169 ENDIF
170
171 dta =half*dt1
172C
173 DO i=1,nel
174 pnew(i) = c1(i)*amu(i) !
175 treps = d1(i)+d2(i)+d3(i)
176 ggt = dt1*gg(i)
177 sig(i,1) = sig(i,1)+two*ggt*d1(i)+lmbd(i)*treps
178 sig(i,2) = sig(i,2)+two*ggt*d2(i)+lmbd(i)*treps
179 sig(i,3) = sig(i,3)+two*ggt*d3(i)+lmbd(i)*treps
180 sig(i,4) = sig(i,4)+ggt*d4(i)
181 sig(i,5) = sig(i,5)+ggt*d5(i)
182 sig(i,6) = sig(i,6)+ggt*d6(i)
183 qold(i) = qnew(i)
184 ENDDO
185 !write(*,'("***** sboltlaw ")')
186 !DO I=1,NEL
187 ! write(*,'(i5,6e10.3)')i,sig(i,1:6)
188 !enddo
189C
190C - To avoid discontinuous element pressure when switching material law
191C cf Law1, Law2, etc <=> Pressure is computed using total formulation P = C1*(rho/rho0-1)
192 IF(ismstr == 1)THEN
193 IF(iresp==0)THEN
194 DO i=1,nel
195 IF(tt >= (bpreld(i,2)-two*dt1)) THEN
196 p = -third*(sig(i,1)+sig(i,2)+sig(i,3))
197 rho(i)=rho0_1*(one+p/c1_1) ! rho will be used in next cycles to compute
198 ! rho_new = rho_old - rho0 dV
199 ! AMU = rho/rho0 - 1 in double precision
200 ENDIF
201 ENDDO
202 ELSE ! Single precision
203 DO i=1,nel
204 IF(tt >= (bpreld(i,2)-two*dt1)) THEN
205 p = -third*(sig(i,1)+sig(i,2)+sig(i,3))
206 rho(i)=rho0_1*(one+p/c1_1)
207 vol0dp(i) = vol0(i)*(one+p/c1_1) ! VOL0DP will be used in next cycle to compute
208 ! AMU = VOL0DP/VOLDP - 1 in single precision, VOLDP=V0 in small strain
209 ENDIF
210 ENDDO
211 END IF
212 ELSEIF(ismstr==2)THEN
213 IF(iresp==0)THEN
214 DO i=1,nel
215 IF(offg(i) > one)THEN
216 IF(tt >= (bpreld(i,2)-two*dt1)) THEN
217 p = -third*(sig(i,1)+sig(i,2)+sig(i,3))
218 rho(i)=rho0_1*(one+p/c1_1) ! rho will be used in next cycles to compute
219 ! rho_new = rho_old - rho0 dV
220 ! AMU = rho/rho0 - 1 in double precision
221 ENDIF
222 ELSE
223 p = -third*(sig(i,1)+sig(i,2)+sig(i,3))
224 vol0(i) = vnew(i)*(one+p/c1_1)
225 rho(i) = rho0_1*vol0(i)/vnew(i) ! cf rho = rho0 * V0 / V at each time step
226 ! amu = rho/rho0 - 1 in double precision
227
228 END IF
229 ENDDO
230 ELSE ! Single precision
231 DO i=1,nel
232 IF(offg(i) > one)THEN
233 IF(tt >= (bpreld(i,2)-two*dt1)) THEN
234 p = -third*(sig(i,1)+sig(i,2)+sig(i,3))
235 rho(i) = rho0_1*(one+p/c1_1)
236 vol0dp(i) = vol0(i)*(one+p/c1_1) ! VOL0DP will be used in next cycle to compute
237 ! AMU = VOL0DP/VOLDP - 1 in single precision, VOLDP=V0 in small strain
238 ENDIF
239 ELSE
240 p = -third*(sig(i,1)+sig(i,2)+sig(i,3))
241 vol0(i) = vnew(i)*(one+p/c1_1)
242 vol0dp(i) = vol0(i) ! cf rho = rho0 * VOL0DP / VOLDP
243 ! amu = vol0dp/voldp - 1 in single precision, voldp = v in large strain
244 rho(i) = rho0_1*vol0(i)/vnew(i)
245 END IF
246 ENDDO
247 END IF
248 ELSEIF(ismstr<=4)THEN
249 IF(iresp==0)THEN
250 DO i=1,nel
251 p = -third*(sig(i,1)+sig(i,2)+sig(i,3))
252 vol0(i)=vnew(i)*(one+p/c1_1)
253 rho(i)= rho0_1*vol0(i)/vnew(i) ! cf rho = rho0 * V0 / V at each time step
254 ! AMU = rho/rho0 - 1 in double precision
255 ENDDO
256 ELSE ! Single precision
257 DO i=1,nel
258 p = -third*(sig(i,1)+sig(i,2)+sig(i,3))
259 vol0(i) = vnew(i)*(one+p/c1_1)
260 rho(i) = rho0_1*vol0(i)/vnew(i)
261 vol0dp(i) = vol0(i) ! cf rho = rho0 * VOL0DP / VOLDP
262 ! AMU = VOL0DP/VOLDP - 1 in single precision, VOLDP = V in large strain
263 ENDDO
264 END IF
265 END IF
266C
267 DO i=1,nel
268 e1=d1(i)*(sold1(i)+sig(i,1))
269 e2=d2(i)*(sold2(i)+sig(i,2))
270 e3=d3(i)*(sold3(i)+sig(i,3))
271 e4=d4(i)*(sold4(i)+sig(i,4))
272 e5=d5(i)*(sold5(i)+sig(i,5))
273 e6=d6(i)*(sold6(i)+sig(i,6))
274 einc= zero !VOL_AVG(I)*(E1+E2+E3+E4+E5+E6)*DTA - HALF*DVOL(I)*(QOLD(I)+QNEW(I))
275 eint(i)=(eint(i)+einc*off(i)) / max(em15,vol0(i))
276 ENDDO
277C
278 RETURN
279 END
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
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
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
sboltlaw
subroutine sboltlaw(pm, sig, mat, d1, d2, d3, d4, d5, d6, nel, rho, bpreld, eint, qold, vol0, stifn, dt2t, neltst, ityptst, offg, geo, pid, mumax, ngl, ssp, dvol, aire, vnew, vd2, deltax, vis, pnew, psh, qnew, ssp_eq, sold1, sold2, sold3, sold4, sold5, sold6, mssa, dmels, conde, amu, vol_avg, dtel, g_dt, off, ipm, rhoref, rhosp, vol0dp, ismstr, jsph, jtur, ity, jthe, jsms, npg, glob_therm, numgeo, igeo)
Definition sboltlaw.F:51