OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps106.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!|| sigeps106 ../engine/source/materials/mat/mat106/sigeps106.F
25!||--- called by ------------------------------------------------------
26!|| mulaw ../engine/source/materials/mat_share/mulaw.F90
27!||--- calls -----------------------------------------------------
28!|| finter ../engine/source/tools/curve/finter.F
29!|| jc ../engine/source/materials/mat/mat106/sigeps106.F
30!||====================================================================
31 SUBROUTINE sigeps106(
32 A NEL , NUPARAM, NUVAR , NFUNC , IFUNC ,
33 B NPF , TF , TIME , TIMESTEP, UPARAM ,
34 C RHO0 , RHO , VOLUME , EINT ,
35 D DEPSXX , DEPSYY , DEPSZZ , DEPSXY , DEPSYZ , DEPSZX ,
36 E EPSXX , EPSYY , EPSZZ , EPSXY , EPSYZ , EPSZX ,
37 F SIGOXX , SIGOYY , SIGOZZ , SIGOXY , SIGOYZ , SIGOZX ,
38 G SIGNXX , SIGNYY , SIGNZZ , SIGNXY , SIGNYZ , SIGNZX ,
39 H SIGVXX , SIGVYY , SIGVZZ , SIGVXY , SIGVYZ , SIGVZX ,
40 I SOUNDSP, VISCMAX, UVAR , OFF , PLA , DPLA ,
41 J EPSD , TEMP , JTHE , JLAG , FHEAT )
42C------------------------------------------------------------------------
43C Radial return for Johnson-Cook without viscous effect
44C R = (A+B*p^n)*(1-Tr^m)
45C--------------------------------------------------
46C Johnson-Cook without viscous effect
47C sigma = (A+B*p^n)*(1-Tr^m) with Tr=(T-Tref)/(Tm-Tref)
48C p - cumulated plastic strain
49C epsm - numerical set max value for plastic strain
50C sigmam - numerical set max value for stress
51C--------------------------------------------------
52C Be careful,
53C
54C EPSXX ... are updates strain values = eps^{n+1}
55C with eps_xx^{n+1} = eps_xx^{n} + deps_xx ...
56
57C-----------------------------------------------
58C I m p l i c i t T y p e s
59C-----------------------------------------------
60#include "implicit_f.inc"
61C----------------------------------------------------------------
62C I n p u t A r g u m e n t s
63C----------------------------------------------------------------
64 INTEGER NEL, NUPARAM, NUVAR
65 INTEGER ,INTENT(IN) :: JTHE
66 INTEGER ,INTENT(IN) :: JLAG
67 my_real
68 . TIME , TIMESTEP , UPARAM(NUPARAM),
69 . RHO(NEL) , RHO0(NEL) , VOLUME(NEL), EINT(NEL),
70 . DEPSXX(NEL), DEPSYY(NEL), DEPSZZ(NEL), DEPSXY(NEL), DEPSYZ(NEL), DEPSZX(NEL),
71 . EPSXX (NEL), EPSYY (NEL), EPSZZ (NEL), EPSXY (NEL), EPSYZ (NEL), EPSZX (NEL),
72 . SIGOXX(NEL), SIGOYY(NEL), SIGOZZ(NEL), SIGOXY(NEL), SIGOYZ(NEL), SIGOZX(NEL),
73 . OFF(NEL) , PLA(NEL) , DPLA(NEL) , EPSD(NEL)
74C----------------------------------------------------------------
75C O u t p u t A r g u m e n t s
76C----------------------------------------------------------------
77 my_real
78 . signxx(nel), signyy(nel), signzz(nel),
79 . signxy(nel), signyz(nel), signzx(nel),
80 . sigvxx(nel), sigvyy(nel), sigvzz(nel),
81 . sigvxy(nel), sigvyz(nel), sigvzx(nel),
82 . soundsp(nel), viscmax(nel)
83C----------------------------------------------------------------
84C I n p u t O u t p u t A r g u m e n t s
85C----------------------------------------------------------------
86 my_real :: uvar(nel,nuvar)
87 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: temp
88 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: fheat
89C----------------------------------------------------------------------
90C V a r i a b l e s f o r f u n c t i o n i n t e r p o l a t i o n
91C----------------------------------------------------------------------
92 INTEGER :: NPF(*), NFUNC, IFUNC(NFUNC)
93 my_real :: TF(*)
94C----------------------------------------------------------------
95C L o c a l V a r i a b l e s
96C----------------------------------------------------------------
97 INTEGER I,J,NMAX,N, ITER, ITER_MAX
98 my_real
99 . e,nu,g,a,b,cm,cn,tm,tref,epsm,sigmam,
100 . mu,mu0,lambda0,mu_ratio,bulk,lambda,
101 . f,tol,interval,e0, nu0,t0,
102 . sxx,syy,szz,sxy,syz,szx,eqvm,trc,ratio,
103 . sigxx,sigyy,sigzz,sigxy,sigyz,sigzx,
104 . exx,eyy,ezz,exy,eyz,ezx,trace_deps,trace_sig, pressure, pressure0,
105 . v0,v1,fun,dfun, scale, cs, error, jc_val, djc_val, incr, alpha,
106 . xx, dxx, norm_fun0, error1, error2, norm_v0,
107 . depsp_xx, depsp_yy, depsp_zz, depsp_xy, depsp_yz, depsp_zx, trceps
108 my_real
109 . finter,dydx
110 LOGICAL :: CONVERGED
111 EXTERNAL FINTER
112C-----------------------------------------------
113C PARAMETERS READING
114C-----------------------------------------------
115 E = uparam(1)
116 nu = uparam(2)
117 a = uparam(3)
118 b = uparam(4)
119 cm = uparam(5)
120 cn = uparam(6)
121 tm = uparam(7)
122 tref = uparam(8)
123 epsm = uparam(11)
124 sigmam = uparam(12)
125 nmax = uparam(13)
126 tol = uparam(14)
127 cs = uparam(15)
128C-----------------------------------------------
129C USER VARIABLES INITIALIZATION
130C-----------------------------------------------
131 IF(timestep == zero)THEN
132 DO i=1,nel
133 uvar(i,1)=zero ! accumulated plastic strain
134 uvar(i,2)=zero ! JC yield stress
135 uvar(i,3)=zero ! Tangent of yield stress
136 uvar(i,4)=zero ! Plastic strain XX
137 uvar(i,5)=zero ! Plastic strain YY
138 uvar(i,6)=zero ! Plastic strain ZZ
139 uvar(i,7)=zero ! Plastic strain XY
140 uvar(i,8)=zero ! Plastic strain YZ
141 uvar(i,9)=zero ! Plastic strain ZX
142C Temperature dependent material properties
143 IF (ifunc(1) > 0) THEN
144 e = uparam(1)*finter(ifunc(1),temp(i),npf,tf,dydx)
145 ENDIF
146 IF (ifunc(3) > 0) THEN
147 nu = uparam(2)*finter(ifunc(3),temp(i),npf,tf,dydx)
148 ENDIF
149 uvar(i,10)=e ! Previous Young's modulus
150 uvar(i,11)=nu ! Previous Poisson's ratio
151 uvar(i,12)=temp(i) ! Previous temperature
152 pla(i) = zero
153 ENDDO
154 ENDIF
155
156 interval=one
157 DO i=1,nel
158 t0 = uvar(i,12) ! previous temperature
159C Temperature dependent material properties
160 IF (ifunc(1) > 0) THEN
161 IF (ifunc(2) == 0 .OR. temp(i) > t0) THEN
162 e = uparam(1)*finter(ifunc(1),temp(i),npf,tf,dydx)
163 ELSEIF (ifunc(2) > 0) THEN
164 e = uparam(1)*finter(ifunc(2),temp(i),npf,tf,dydx)
165 ENDIF
166 ENDIF
167 IF (ifunc(3) > 0) THEN
168 nu = uparam(2)*finter(ifunc(3),temp(i),npf,tf,dydx)
169 ENDIF
170C Important, values newly added element should be initialized
171 IF (uvar(i,10) == zero) THEN
172 uvar(i,10) = e
173 ENDIF
174 IF (uvar(i,11) == zero) THEN
175 uvar(i,11) = nu
176 ENDIF
177
178C E and NU at previous time
179 nu0 = uvar(i, 11)
180 e0 = uvar(i, 10)
181
182 mu = half * e / (one + nu)
183 mu0 = half * e0 / (one + nu0)
184 mu_ratio = mu / mu0
185 lambda = e * nu / (one + nu) / (one - two * nu)
186 lambda0 = e0 * nu0 / (one + nu0) / (one - two * nu0)
187 bulk = third * e / (one - two * nu)
188
189C Sound speed = sqrt((lambda+2*mu)/density)
190 soundsp(i) = sqrt((two*mu+lambda)/rho(i))
191 viscmax(i) = zero
192
193C Deviatoric stress at previous time step
194 trc = sigoxx(i) + sigoyy(i) + sigozz(i)
195 pressure0 = - third * trc
196 sxx = sigoxx(i) + pressure0
197 syy = sigoyy(i) + pressure0
198 szz = sigozz(i) + pressure0
199 sxy = sigoxy(i)
200 syz = sigoyz(i)
201 szx = sigozx(i)
202
203C Deviatoric strain increment
204 trceps = depsxx(i) + depsyy(i) + depszz(i)
205 exx = depsxx(i) - third * trceps
206 eyy = depsyy(i) - third * trceps
207 ezz = depszz(i) - third * trceps
208 exy = depsxy(i)
209 eyz = depsyz(i)
210 ezx = depszx(i)
211
212C Deviatoric elastic prediction at current time step
213 sxx = (sxx * mu_ratio + two * mu * exx)
214 syy = (syy * mu_ratio + two * mu * eyy)
215 szz = (szz * mu_ratio + two * mu * ezz)
216 sxy = (sxy * mu_ratio + mu * exy)
217 syz = (syz * mu_ratio + mu * eyz)
218 szx = (szx * mu_ratio + mu * ezx)
219C Trace of the updated strain tensor
220 trc = epsxx(i) + epsyy(i) + epszz(i)
221C Pressure at time n+1 does not depend on plastic strain ...
222C a direct computation leads to:
223 pressure = - bulk * trc
224C=======================================================================
225C II - JC yield stress
226C=======================================================================
227 CALL jc(uvar(i,1),temp(i),a,b,cm,cn,tref,tm,epsm,sigmam,uvar(i,2),uvar(i,3))
228 eqvm = sqrt(three_half*(sxx**2+syy**2+szz**2+two*(sxy**2+syz**2+szx**2)))
229 f = eqvm - uvar(i,2)
230
231 uvar(i, 13) = zero
232 uvar(i, 14) = zero
233
234 IF (f <= zero) THEN
235C=======================================================================
236C Case 1: pure elastic
237C=======================================================================
238 dpla(i) = zero
239 ratio = one
240 ELSEIF(temp(i) >= tm) THEN
241C=======================================================================
242C Case 2: pure hydrodynamic
243C=======================================================================
244 dpla(i) = zero
245 ratio = zero
246 ELSE
247C=======================================================================
248C Case 3: elastoplastic
249C Newton-Raphson to solve:
250C eqVM - 3*mu*dp - R(p0+dp)=0
251C to get dp
252C=======================================================================
253C Initial condition for Newton method
254 v0 = eqvm / (three * mu)
255 CALL jc(uvar(i,1)+v0,temp(i),a,b,cm,cn,tref,tm,epsm,sigmam,jc_val,djc_val)
256 fun = eqvm - three * mu * v0 - jc_val
257 norm_fun0 = abs(fun)
258 norm_v0 = abs(v0)
259 converged = .false.
260 iter = 0
261 error2 = zero
262 iter_max = nmax
263 DO WHILE (.NOT. converged .AND. iter <= iter_max)
264 CALL jc(uvar(i,1)+v0,temp(i),a,b,cm,cn,tref,tm,epsm,sigmam,jc_val,djc_val)
265 fun = eqvm - three * mu * v0 - jc_val
266 error1 = abs(fun)
267 converged = error1 < tol * norm_fun0
268 IF (.NOT. converged) THEN
269 dfun = -three * mu - djc_val
270 incr = - fun / dfun
271 alpha = one
272 IF (incr < zero) THEN
273 alpha =min(one , - nine / ten * v0 / incr)
274 ENDIF
275 v0 = v0 + alpha * incr
276 error2 = abs(alpha * incr)
277 converged = error2 < tol * abs(v0)
278 ENDIF
279 iter = iter + 1
280 ENDDO
281
282 uvar(i, 13) = iter
283 uvar(i, 14) = error2
284
285 dpla(i) = v0
286 uvar(i,1) = uvar(i,1) + v0
287 CALL jc(uvar(i,1),temp(i),a,b,cm,cn,tref,tm,epsm,sigmam,uvar(i,2),uvar(i,3))
288 scale=three_half*v0/eqvm
289 uvar(i,4) = uvar(i,4) + scale * sxx
290 uvar(i,5) = uvar(i,5) + scale * syy
291 uvar(i,6) = uvar(i,6) + scale * szz
292 uvar(i,7) = uvar(i,7) + scale * sxy
293 uvar(i,8) = uvar(i,8) + scale * syz
294 uvar(i,9) = uvar(i,9) + scale * szx
295
296 uvar(i,2)=uvar(i,2)
297 uvar(i,3)=uvar(i,3)
298 ratio = (one - three * mu * v0 / eqvm)
299!
300 IF (jthe /= 0 .AND. jlag /= 0) THEN
301 fheat(i) = fheat(i) + uvar(i,2)*dpla(i)*volume(i)
302 ELSE IF (cs > zero) THEN
303 temp(i) = temp(i) + uvar(i,2)*dpla(i)/cs
304 ENDIF
305 ENDIF
306
307 uvar(i,10) = e ! register as previous Young's modulus
308 uvar(i,11) = nu ! register as previous Poisson's ratio
309 uvar(i,12) = temp(i) ! register as previous temperature
310 signxx(i) = sxx * ratio - pressure
311 signyy(i) = syy * ratio - pressure
312 signzz(i) = szz * ratio - pressure
313 signxy(i) = sxy * ratio
314 signyz(i) = syz * ratio
315 signzx(i) = szx * ratio
316 ENDDO
317C
318 pla(1:nel) = uvar(1:nel,1)
319
320 DO i=1,nel
321 viscmax(i)= zero
322 sigvxx(i) = zero
323 sigvyy(i) = zero
324 sigvzz(i) = zero
325 sigvxy(i) = zero
326 sigvyz(i) = zero
327 sigvzx(i) = zero
328 ENDDO
329
330 RETURN
331 END
332
333!||====================================================================
334!|| jc ../engine/source/materials/mat/mat106/sigeps106.F
335!||--- called by ------------------------------------------------------
336!|| sigeps106 ../engine/source/materials/mat/mat106/sigeps106.F
337!||====================================================================
338 SUBROUTINE jc(P,T,A,B,CM,CN,TREF,TM,EPSM,SIGMAM,JC_YIELD,TAN_JC)
339C--------------------------------------------------
340C Johnson-Cook and tangent without viscous effect
341C--------------------------------------------------
342C-----------------------------------------------
343C I m p l i c i t T y p e s
344C-----------------------------------------------
345#include "implicit_f.inc"
346C-----------------------------------------------
347C D u m m y A r g u m e n t s
348C-----------------------------------------------
349 my_real, INTENT(IN) :: p,t,a,b,cm,cn
350 my_real, INTENT(IN) :: tref,tm,epsm,sigmam
351 my_real, INTENT(OUT) :: jc_yield, tan_jc
352 my_real :: tr
353C---------------------------------------------------
354 tr=max(t,tref)
355 tr=min(tr,tm)
356 tr=(tr-tref)/(tm-tref)
357 IF (tm == zero) THEN
358 tr = zero
359 ENDIF
360 IF(p<=zero) THEN
361 jc_yield=a*(one-tr**cm)
362 tan_jc=zero
363 ELSEIF(p>epsm) THEN
364 jc_yield=(a+b*epsm**cn)*(one-tr**cm)
365 tan_jc=zero
366 ELSE
367 jc_yield=(a+b*p**cn)*(one-tr**cm)
368 tan_jc=b*cn*p**(cn-one)*(one-tr**cm)
369 ENDIF
370 jc_yield=min(sigmam,jc_yield)
371 RETURN
372 END SUBROUTINE jc
my_real
#define my_real
Definition cppsort.cpp:32
alpha
#define alpha
Definition eval.h:35
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
jc
subroutine jc(p, t, a, b, cm, cn, tref, tm, epsm, sigmam, jc_yield, tan_jc)
Definition sigeps106.F:339
sigeps106
subroutine sigeps106(nel, nuparam, nuvar, nfunc, ifunc, npf, tf, time, timestep, uparam, rho0, rho, volume, eint, depsxx, depsyy, depszz, depsxy, depsyz, depszx, epsxx, epsyy, epszz, epsxy, epsyz, epszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, sigvxx, sigvyy, sigvzz, sigvxy, sigvyz, sigvzx, soundsp, viscmax, uvar, off, pla, dpla, epsd, temp, jthe, jlag, fheat)
Definition sigeps106.F:42