OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps120_vm.F File Reference
#include "implicit_f.inc"
#include "units_c.inc"
#include "comlock.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine sigeps120_vm (nel, ngl, nuparam, nuvar, time, timestep, uparam, uvar, off, pla, epsd, soundsp, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, depsxx, depsyy, depszz, depsxy, depsyz, depszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, inloc, dplanl, dmg, dmg_scale)

Function/Subroutine Documentation

◆ sigeps120_vm()

subroutine sigeps120_vm ( integer nel,
integer, dimension(nel), intent(in) ngl,
integer nuparam,
integer nuvar,
time,
timestep,
intent(in) uparam,
intent(inout) uvar,
intent(inout) off,
intent(inout) pla,
intent(out) epsd,
intent(out) soundsp,
intent(in) epspxx,
intent(in) epspyy,
intent(in) epspzz,
intent(in) epspxy,
intent(in) epspyz,
intent(in) epspzx,
intent(in) depsxx,
intent(in) depsyy,
intent(in) depszz,
intent(in) depsxy,
intent(in) depsyz,
intent(in) depszx,
intent(in) sigoxx,
intent(in) sigoyy,
intent(in) sigozz,
intent(in) sigoxy,
intent(in) sigoyz,
intent(in) sigozx,
intent(out) signxx,
intent(out) signyy,
intent(out) signzz,
intent(out) signxy,
intent(out) signyz,
intent(out) signzx,
integer inloc,
intent(inout) dplanl,
intent(inout) dmg,
intent(inout) dmg_scale )

Definition at line 28 of file sigeps120_vm.F.

36C-----------------------------------------------
37C I m p l i c i t T y p e s
38C-----------------------------------------------
39#include "implicit_f.inc"
40C-----------------------------------------------
41C C o m m o n B l o c k s
42C-----------------------------------------------
43#include "units_c.inc"
44#include "comlock.inc"
45C----------------------------------------------------------------
46C D u m m y A R G U M E N T S
47C----------------------------------------------------------------
48 INTEGER :: NEL,NUPARAM,NUVAR,INLOC
49 my_real :: time,timestep
50 INTEGER ,DIMENSION(NEL) ,INTENT(IN) :: NGL
51 my_real ,DIMENSION(NUPARAM),INTENT(IN) :: uparam
52 my_real ,DIMENSION(NEL) ,INTENT(IN) ::
53 . epspxx,epspyy,epspzz,epspxy,epspyz,epspzx,
54 . depsxx,depsyy,depszz,depsxy,depsyz,depszx,
55 . sigoxx,sigoyy,sigozz,sigoxy,sigoyz,sigozx
56c
57 my_real ,DIMENSION(NEL) ,INTENT(OUT) :: epsd,soundsp
58 my_real ,DIMENSION(NEL) ,INTENT(OUT) ::
59 . signxx,signyy,signzz,signxy,signyz,signzx
60 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: off,pla,dplanl,dmg,dmg_scale
61 my_real ,DIMENSION(NEL,NUVAR) ,INTENT(INOUT) :: uvar
62C----------------------------------------------------------------
63C L O C A L V A R I B L E S
64C----------------------------------------------------------------
65 INTEGER :: I,II,ITER,NITER,NINDX,NINDF,ITRX,IDAM
66 INTEGER ,DIMENSION(NEL) :: INDX,INDF
67c material parameters
68 my_real :: young,nu,ssp,g,g2,bulk,yld0,qq,beta,hh,c11,c12,
69 . a1f,a2f,as,d1c,d2c,d1f,d2f,d_trx,d_jc,dm,expn
70c Johnson & Cook rate-dependency
71 my_real cc,epsdref,epsdmax,frate,fcut,alpha,alphai
72c Local variables
73 my_real facr,triax,rvoce,fvm,i1p,ip,jp,tp,fact,facd,dfac,dtinv,
74 . dlam,ddep,dfdlam,dyld_dpla,dpla_dlam,dphi_dlam,
75 . gamc,gamf,np_xx,np_yy,np_zz,np_xy,np_yz,np_zx,
76 . nf_xx,nf_yy,nf_zz,nf_xy,nf_yz,nf_zx,dpxx,dpyy,dpzz,dpxy,dpyz,dpzx
77 my_real ,DIMENSION(NEL) :: i1,j2,yld,phi,dam,jcrate,
78 . sxx,syy,szz,sxy,syz,szx,dpla,jcr_log,gamma,dgamm,pla_nl,dpdt_nl
79c--------------------------------
80c Material parameters
81c--------------------------------
82c UPARAM(1) = Young modulus E
83c UPARAM(2) = Poisson's ratio nu
84c UPARAM(3) = shear modulus G
85c UPARAM(4) = bulk modulus K
86c UPARAM(5) = Yld0: instant yield stress
87c UPARAM(6) = Q : expontial term coefficient in the hardening law
88c UPARAM(7) = BETA : exponent of the exponential term
89c UPARAM(8) = P : multiplier of the linear term in the hardening law
90c UPARAM(9) = A1F : parameter of the yield function
91c UPARAM(10) = A2F : parameter of the yield function
92c UPARAM(11) = A1H : distortionnal yield hardening coefficiant
93c UPARAM(12) = A2H : distortionnal yield hardening coefficiant
94c UPARAM(13) = AS : parameter of the potential function
95c UPARAM(14) = D1C: first damage strain parameter in initial damage threshold
96c UPARAM(15) = D2C: second damage strain parameter in initial damage threshold
97c UPARAM(16) = D1F: first damage strain parameter in final damage threshold
98c UPARAM(17) = D2F: second damage strain parameter in final damage threshold
99c UPARAM(18) = D_TRX : triaxiality factor in damage formula
100c UPARAM(19) = D_JC: JC strain rate coefficient in damage formula
101c UPARAM(20) = EXPN exponent in damage evolution
102c UPARAM(21) = CC : Johnson-Cook strain rate-dependency coefficient
103c UPARAM(22) = EPSDREF quasi-static reference strain rate
104c UPARAM(23) = EPSDMAX maximal reference strain rate
105c UPARAM(24) = FCUT : cut frequency for strain rate filtering
106c UPARAM(25) = IFORM = 1: Yield formulation flag => Drucker-Prager in tension
107c IFORM = 2: Yield formulation flag => Von Mises in tension
108c UPARAM(26) = ITRX = 1 : pressure dependent for all T values
109c ITRX = 2 : no pressure dependency when T < 0
110c UPARAM(27) = IDAM = 1 : damage model without turning point
111c IDAM = 2 : damage model with turning point
112c UPARAM(28) = SSP : sound speed
113c UPARAM(29) = Table Id
114c UPARAM(30) = Xscale for yld function
115c UPARAM(31) = Yscale for yld function
116c UPARAM(32) = Elastic modulus C11
117c UPARAM(33) = Elastic modulus C12
118c--------------------------------
119c UVAR(1) : ARCL = PLA / (1-DAM) or non local plastic strain
120c UVAR(2) : DAM
121c UVAR(3) : TOTAL STRAIN RATE
122c UVAR(4) : plastic function residual (for output/convergence check)
123c UVAR(5) : non local plastic strain
124C=======================================================================
125 young = uparam(1)
126 nu = uparam(2)
127 g = uparam(3)
128 bulk = uparam(4)
129c Parameters of yield function/hardening law and flow potential
130 yld0 = uparam(5)
131 qq = uparam(6)
132 beta = uparam(7)
133 hh = uparam(8)
134 a1f = uparam(9)
135 a2f = uparam(10)
136c A1H = UPARAM(11)
137c A2H = UPARAM(12)
138 as = uparam(13)
139c Johnson-Cook failure parameters
140 d1c = uparam(14)
141 d2c = uparam(15)
142 d1f = uparam(16)
143 d2f = uparam(17)
144 d_trx = uparam(18)
145 d_jc = uparam(19)
146 expn = uparam(20)
147c Johnson-Cook strain rate-dependency and distortional hardening
148 cc = uparam(21)
149 epsdref = uparam(22)
150 epsdmax = uparam(23)
151 fcut = uparam(24)
152c
153 itrx = nint(uparam(26))
154 idam = nint(uparam(27))
155 ssp = uparam(28)
156 c11 = uparam(32)
157 c12 = uparam(33)
158c
159 g2 = g * two
160 alpha = 0.005
161 alphai = one - alpha
162c
163 IF (inloc > 0) THEN
164 dtinv = one / max(timestep, em20)
165 DO i = 1,nel
166 pla_nl(i) = uvar(i,5) + max(dplanl(i),zero)
167 uvar(i,5) = pla_nl(i)
168 dpdt_nl(i) = min(max(dplanl(i),zero)*dtinv ,epsdmax)
169 ENDDO
170 ENDIF
171c------------------------------------------------------------------
172c Computation of the nominal trial stress tensor (non damaged)
173c------------------------------------------------------------------
174 dam(1:nel) = dmg(1:nel)
175 dpla(1:nel) = zero
176c
177 DO i = 1,nel
178 IF (off(i) < 0.001) off(i) = zero
179 IF (off(i) < one) off(i) = off(i)*four_over_5
180 IF (off(i) == one) THEN
181 signxx(i) = sigoxx(i) + c11*depsxx(i) + c12*(depsyy(i) + depszz(i))
182 signyy(i) = sigoyy(i) + c11*depsyy(i) + c12*(depsxx(i) + depszz(i))
183 signzz(i) = sigozz(i) + c11*depszz(i) + c12*(depsxx(i) + depsyy(i))
184 signxy(i) = sigoxy(i) + depsxy(i)*g
185 signyz(i) = sigoyz(i) + depsyz(i)*g
186 signzx(i) = sigozx(i) + depszx(i)*g
187 i1(i) = signxx(i) + signyy(i) + signzz(i)
188c
189 sxx(i) = signxx(i) - i1(i) * third
190 syy(i) = signyy(i) - i1(i) * third
191 szz(i) = signzz(i) - i1(i) * third
192 sxy(i) = signxy(i)
193 syz(i) = signyz(i)
194 szx(i) = signzx(i)
195 j2(i) = (sxx(i)**2 + syy(i)**2 + szz(i)**2)*half
196 . + sxy(i)**2 + syz(i)**2 + szx(i)**2
197 END IF
198 ENDDO
199c Johnson & Cook rate dependency (total strain rate)
200 jcr_log(1:nel) = zero
201 jcrate(1:nel) = one
202 IF (epsdref > zero) THEN
203 IF (inloc == 1) THEN ! non local plastic strain rate
204 jcr_log(1:nel) = log(dpdt_nl(1:nel) / epsdref)
205 jcrate(1:nel) = one + cc * jcr_log(1:nel)
206 ELSE ! total strain rate
207 DO i = 1,nel
208 IF (off(i) == one) THEN
209 epsd(i) = (epspxx(i)**2 + epspyy(i)**2 + epspzz(i)**2 )*two
210 . + epspxy(i)**2 + epspyz(i)**2 + epspzx(i)**2
211 epsd(i) = min(sqrt(epsd(i)) ,epsdmax)
212 epsd(i) = alpha*epsd(i) + alphai*uvar(i,3)
213 uvar(i,3) = epsd(i)
214 IF (epsd(i) > epsdref) THEN
215 jcr_log(i) = log(epsd(i) / epsdref)
216 jcrate(i) = one + cc * jcr_log(i)
217 END IF
218 END IF
219 ENDDO
220 END IF
221 END IF
222c------------------------------------------------------------------
223c Actual yield value, yield criterion and plastic flow function
224c------------------------------------------------------------------
225 DO i = 1,nel
226 rvoce = qq*(one - exp(-beta*pla(i))) + hh*pla(i)
227 yld(i) = (yld0 + rvoce) * jcrate(i)
228 fvm = max(zero, i1(i) + half*sqr3*a1f*yld0 / a2f )
229 phi(i) = j2(i) + third*a2f * fvm**2
230 . - (yld(i)**2 + fourth*(a1f*yld0)**2/a2f)
231 ENDDO
232c------------------------------------------------------------------
233c Plasticity
234c------------------------------------------------------------------
235 nindx = 0
236 nindf = 0
237 DO i = 1,nel
238 IF (phi(i) > zero .and. off(i) == one) THEN
239 nindx = nindx + 1
240 indx(nindx) = i
241 END IF
242 ENDDO
243c
244 DO ii = 1,nindx
245 i = indx(ii)
246 dpla(i) = zero
247 i1p = max(zero, i1(i))
248 fvm = max(zero, i1(i) + half*sqr3*a1f*yld0 / a2f)
249 niter = 4
250c
251 DO iter = 1,niter
252
253 ! normal to non associated plastic flow surface = d_psi/d_sig
254 jp = third * as * i1p
255 tp = two * jp
256c
257 np_xx = sxx(i) + tp
258 np_yy = syy(i) + tp
259 np_zz = szz(i) + tp
260 np_xy = sxy(i)
261 np_yz = syz(i)
262 np_zx = szx(i)
263c
264 ! normal to plastic yield surface = d_phi/d_sig
265 ip = two_third * a2f * fvm
266 nf_xx = sxx(i) + ip
267 nf_yy = syy(i) + ip
268 nf_zz = szz(i) + ip
269 nf_xy = sxy(i)
270 nf_yz = syz(i)
271 nf_zx = szx(i)
272c
273 ! derivative of yld criterion : dphi/dlam = d_phi/d_sig * dsig/dlam
274 dfdlam =-nf_xx * (c11*np_xx + c12 * (np_yy + np_zz))
275 . - nf_yy * (c11*np_yy + c12 * (np_xx + np_zz))
276 . - nf_zz * (c11*np_zz + c12 * (np_xx + np_yy))
277 . -(nf_xy*np_xy + nf_yz*np_yz + nf_zx*np_zx)*two*g2
278c
279 dyld_dpla = (hh + beta*qq*exp(-beta*pla(i))) * jcrate(i)
280 dpla_dlam = two*sqrt((j2(i) + tp*as*i1(i)))
281 dphi_dlam = dfdlam - two*yld(i)*dyld_dpla*dpla_dlam
282c----------------------------------------------------------------
283 dlam = -phi(i) / dphi_dlam
284c----------------------------------------------------------------
285 ! Plastic strains tensor update
286 dpxx = dlam*np_xx
287 dpyy = dlam*np_yy
288 dpzz = dlam*np_zz
289 dpxy = dlam*np_xy
290 dpyz = dlam*np_yz
291 dpzx = dlam*np_zx
292 ! Elasto-plastic stresses update
293 signxx(i) = signxx(i) - (c11*dpxx + c12*(dpyy + dpzz))
294 signyy(i) = signyy(i) - (c11*dpyy + c12*(dpxx + dpzz))
295 signzz(i) = signzz(i) - (c11*dpzz + c12*(dpxx + dpyy))
296 signxy(i) = signxy(i) - g2*dpxy
297 signyz(i) = signyz(i) - g2*dpyz
298 signzx(i) = signzx(i) - g2*dpzx
299c
300 ! Plastic strain increments update
301 ddep = dlam * dpla_dlam
302 dpla(i)= max(zero, dpla(i) + ddep)
303 pla(i) = pla(i) + ddep
304c----------------------
305c criterion update
306c----------------------
307 i1(i) = signxx(i) + signyy(i) + signzz(i)
308 i1p = max(zero, i1(i))
309 sxx(i) = signxx(i) - i1(i)*third
310 syy(i) = signyy(i) - i1(i)*third
311 szz(i) = signzz(i) - i1(i)*third
312 sxy(i) = signxy(i)
313 syz(i) = signyz(i)
314 szx(i) = signzx(i)
315 j2(i) =(sxx(i)**2 + syy(i)**2 + szz(i)**2 ) * half
316 . + sxy(i)**2 + syz(i)**2 + szx(i)**2
317 rvoce = qq*(one - exp(-beta*pla(i))) + hh*pla(i)
318 yld(i) = (yld0 + rvoce) * jcrate(i)
319 fvm = max(zero, i1(i) + half*sqr3*a1f*yld0 / a2f)
320 phi(i) = j2(i) + third*a2f * fvm**2
321 . - (yld(i)**2 + fourth*(a1f*yld0)**2/a2f)
322 END DO ! Newton iterations
323c
324 uvar(i,4) = phi(i) / yld0**2
325 ENDDO ! II = 1,NINDX
326c---------------------------------------------------------------------
327c Update damage
328c---------------------------------------------------------------------
329 IF (idam > 0) THEN
330 IF (idam == 1) THEN
331 IF (inloc == 1) THEN
332 dgamm(1:nel) = dplanl(1:nel)
333 gamma(1:nel) = pla_nl(1:nel)
334 ELSE
335 dgamm(1:nel) = dpla(1:nel)
336 gamma(1:nel) = pla(1:nel)
337 END IF
338 ELSE
339 IF (inloc == 1) THEN
340 dgamm(1:nel) = dplanl(1:nel) * dmg_scale(1:nel)
341 uvar(1:nel,1) = uvar(1:nel,1) + dgamm(1:nel)
342 gamma(1:nel) = uvar(1:nel,1)
343 ELSE
344 dgamm(1:nel) = dpla(1:nel) * dmg_scale(1:nel)
345 uvar(1:nel,1) = uvar(1:nel,1) + dgamm(1:nel)
346 gamma(1:nel) = uvar(1:nel,1)
347 END IF
348 END IF
349c
350 DO i = 1,nel
351c
352 triax = zero
353 fact = one
354 facr = one + d_jc * jcr_log(i) ! total strain rate factor on damage
355 IF ( j2(i)/= zero)triax = third * i1(i) / sqrt(three*j2(i))
356 IF (itrx == 1 ) THEN
357 fact = exp(-d_trx*triax)
358 ELSE
359 IF (triax > zero )fact = exp(-d_trx*triax)
360 ENDIF
361 gamc = (d1c + d2c * fact) * facr
362 gamf = (d1f + d2f * fact) * facr
363 IF (gamma(i) > gamc) THEN
364 IF (dam(i) == zero) THEN
365 nindf = nindf + 1
366 indf(nindf) = i
367 END IF
368 IF (expn == one) THEN
369 dam(i) = dam(i) + dgamm(i) / (gamf - gamc)
370 ELSE
371 dfac = (gamma(i) - gamc) / (gamf - gamc)
372 dam(i) = dam(i) + expn * dfac**(expn-one) * dgamm(i) / (gamf - gamc)
373 END IF
374 IF (dam(i) >= one) THEN
375 dam(i) = one
376 off(i) = four_over_5
377 END IF
378 dmg(i) = dam(i)
379 END IF ! GAMMA > GAMC
380 ENDDO
381c Update damaged stress
382 DO i = 1, nel
383 dmg_scale(i) = one - dam(i)
384 END DO
385 END IF
386c-------------------------
387 soundsp(1:nel) = ssp
388c-------------------------
389 IF (nindf > 0) THEN
390 DO ii=1,nindf
391#include "lockon.inc"
392 WRITE(iout, 1000) ngl(indf(ii))
393 WRITE(istdo,1100) ngl(indf(ii)),time
394#include "lockoff.inc"
395 ENDDO
396 END IF
397c-----------------------------------------------------------------
398 1000 FORMAT(1x,'START DAMAGE IN SOLID ELEMENT NUMBER ',i10)
399 1100 FORMAT(1x,'START DAMAGE IN SOLID ELEMENT NUMBER ',i10,1x,' AT TIME :',g11.4)
400c-----------------------------------------------------------------
401 RETURN
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