OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps120_dp.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_dp (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_dp()

subroutine sigeps120_dp ( 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_dp.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,a1h,a2h,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,fdp,i1p,ip,jp,tp,fact,facd,dfac,dtinv,
74 . dlam,ddep,dfdlam,dfdpla,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,a1,a2,yld,phi,dam,fdam,jcrate,
78 . sxx,syy,szz,sxy,syz,szx,dpla,jcr_log,gamma,gammav,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)
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)
136 a1h = uparam(11)
137 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
176
177c
178 DO i = 1,nel
179 IF (off(i) < 0.001) off(i) = zero
180 IF (off(i) < one) off(i) = off(i)*four_over_5
181 IF (off(i) == one) THEN
182 signxx(i) = sigoxx(i) + c11*depsxx(i) + c12*(depsyy(i) + depszz(i))
183 signyy(i) = sigoyy(i) + c11*depsyy(i) + c12*(depsxx(i) + depszz(i))
184 signzz(i) = sigozz(i) + c11*depszz(i) + c12*(depsxx(i) + depsyy(i))
185 signxy(i) = sigoxy(i) + depsxy(i)*g
186 signyz(i) = sigoyz(i) + depsyz(i)*g
187 signzx(i) = sigozx(i) + depszx(i)*g
188 i1(i) = signxx(i) + signyy(i) + signzz(i)
189c
190 sxx(i) = signxx(i) - i1(i) * third
191 syy(i) = signyy(i) - i1(i) * third
192 szz(i) = signzz(i) - i1(i) * third
193 sxy(i) = signxy(i)
194 syz(i) = signyz(i)
195 szx(i) = signzx(i)
196 j2(i) = (sxx(i)**2 + syy(i)**2 + szz(i)**2)*half
197 . + sxy(i)**2 + syz(i)**2 + szx(i)**2
198 END IF
199 ENDDO
200c Johnson & Cook rate dependency
201 jcr_log(1:nel) = zero
202 jcrate(1:nel) = one
203 IF (epsdref > zero) THEN
204 IF (inloc == 1) THEN ! non local plastic strain rate
205 jcr_log(1:nel) = log(dpdt_nl(1:nel) / epsdref)
206 jcrate(1:nel) = one + cc * jcr_log(1:nel)
207 ELSE ! total strain rate
208 DO i = 1,nel
209 IF (off(i) == one) THEN
210 epsd(i) = (epspxx(i)**2 + epspyy(i)**2 + epspzz(i)**2 )*two
211 . + epspxy(i)**2 + epspyz(i)**2 + epspzx(i)**2
212 epsd(i) = min(sqrt(epsd(i)) ,epsdmax)
213 epsd(i) = alpha*epsd(i) + alphai*uvar(i,3)
214 uvar(i,3) = epsd(i)
215 IF (epsd(i) > epsdref) THEN
216 jcr_log(i) = log(epsd(i) / epsdref)
217 jcrate(i) = one + cc * jcr_log(i)
218 END IF
219 END IF
220 ENDDO
221 END IF
222 END IF
223c------------------------------------------------------------------
224c Actual yield value, yield criterion and plastic flow function
225c------------------------------------------------------------------
226 DO i = 1,nel
227 rvoce = qq*(one - exp(-beta*pla(i))) + hh*pla(i)
228 yld(i)= (yld0 + rvoce) * jcrate(i)
229 a1(i) = a1f + a1h * pla(i)
230 a2(i) = max(zero, a2f + a2h * pla(i))
231 fdp = third*(sqr3*yld0*a1(i)*i1(i) + a2(i)*max(zero,i1(i))**2)
232 phi(i)= j2(i) + fdp - yld(i)**2
233 ENDDO
234c------------------------------------------------------------------
235c Plasticity
236c------------------------------------------------------------------
237 nindx = 0
238 nindf = 0
239 DO i = 1,nel
240 IF (phi(i) > zero .and. off(i) == one) THEN
241 nindx = nindx + 1
242 indx(nindx) = i
243 END IF
244 ENDDO
245c
246 DO ii = 1,nindx
247 i = indx(ii)
248 dpla(i) = zero
249 i1p = max(zero,i1(i))
250 niter = 4
251c
252 DO iter = 1,niter
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 = third*(sqr3*yld0*a1(i) + two*a2(i)*i1p)
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 dfdpla = third*(sqr3*a1h*yld0*i1(i) + a2h*i1p**2)
281
282 dpla_dlam = two*sqrt((j2(i) + tp*as*i1(i)))
283c
284 dphi_dlam = dfdlam + (dfdpla - two*yld(i)*dyld_dpla)*dpla_dlam
285c----------------------------------------------------------------
286 dlam = -phi(i) / dphi_dlam
287c----------------------------------------------------------------
288 ! Plastic strains tensor update
289 dpxx = dlam*np_xx
290 dpyy = dlam*np_yy
291 dpzz = dlam*np_zz
292 dpxy = dlam*np_xy
293 dpyz = dlam*np_yz
294 dpzx = dlam*np_zx
295 ! Elasto-plastic stresses update
296 signxx(i) = signxx(i) - (c11*dpxx + c12*(dpyy + dpzz))
297 signyy(i) = signyy(i) - (c11*dpyy + c12*(dpxx + dpzz))
298 signzz(i) = signzz(i) - (c11*dpzz + c12*(dpxx + dpyy))
299 signxy(i) = signxy(i) - g2*dpxy
300 signyz(i) = signyz(i) - g2*dpyz
301 signzx(i) = signzx(i) - g2*dpzx
302c
303 ! Plastic strain increments update
304 ddep = dlam * dpla_dlam
305 dpla(i)= max(zero, dpla(i) + ddep)
306 pla(i) = pla(i) + ddep
307c----------------------
308c criterion update
309c----------------------
310 i1(i) = signxx(i) + signyy(i) + signzz(i)
311 i1p = max(zero, i1(i))
312 sxx(i) = signxx(i) - i1(i)*third
313 syy(i) = signyy(i) - i1(i)*third
314 szz(i) = signzz(i) - i1(i)*third
315 sxy(i) = signxy(i)
316 syz(i) = signyz(i)
317 szx(i) = signzx(i)
318 j2(i) =(sxx(i)**2 + syy(i)**2 + szz(i)**2 ) * half
319 . + sxy(i)**2 + syz(i)**2 + szx(i)**2
320 rvoce = qq*(one - exp(-beta*pla(i))) + hh*pla(i)
321 yld(i) = (yld0 + rvoce) * jcrate(i)
322 a1(i) = a1f + a1h * pla(i)
323 a2(i) = max(zero, a2f + a2h * pla(i))
324 fdp = third*(sqr3*yld0*a1(i)*i1(i) + a2(i)*i1p**2)
325 phi(i) = j2(i) + fdp - yld(i)**2
326 END DO ! Newton iterations
327c
328 !GAMMAV(I) = SQRT(TWO* ( DEPXX(I)**2 + DEPYY(I)**2 + DEPZZ(I)**2))
329 uvar(i,4) = phi(i) / yld0**2
330 ENDDO ! II = 1,NINDX
331c---------------------------------------------------------------------
332c Update damage
333c---------------------------------------------------------------------
334 IF (idam > 0) THEN
335 IF (idam == 1) THEN
336 IF (inloc == 1) THEN
337 dgamm(1:nel) = dplanl(1:nel)
338 gamma(1:nel) = pla_nl(1:nel)
339 ELSE
340 dgamm(1:nel) = dpla(1:nel)
341 gamma(1:nel) = pla(1:nel)
342 END IF
343 ELSE
344 IF (inloc == 1) THEN
345 dgamm(1:nel) = dplanl(1:nel) * dmg_scale(1:nel)
346 uvar(1:nel,1) = uvar(1:nel,1) + dgamm(1:nel)
347 gamma(1:nel) = uvar(1:nel,1)
348 ELSE
349 dgamm(1:nel) = dpla(1:nel) * dmg_scale(1:nel)
350 uvar(1:nel,1) = uvar(1:nel,1) + dgamm(1:nel)
351 gamma(1:nel) = uvar(1:nel,1)
352 END IF
353 END IF
354c
355 DO i = 1,nel ! loop over all element for non local
356
357 triax = zero
358 fact = one
359 facr = one + d_jc * jcr_log(i) ! total strain rate factor on damage
360 IF ( j2(i)/= zero)triax = third * i1(i) / sqrt(three*j2(i))
361 IF (itrx == 1 ) THEN
362 fact = exp(-d_trx*triax)
363 ELSE
364 IF (triax > zero )fact = exp(-d_trx*triax)
365 ENDIF
366 gamc = (d1c + d2c * fact) * facr
367 gamf = (d1f + d2f * fact) * facr
368 IF (gamma(i) > gamc) THEN
369 IF (dam(i) == zero) THEN
370 nindf = nindf + 1
371 indf(nindf) = i
372 END IF
373 IF (expn == one) THEN
374 dam(i) = dam(i) + dgamm(i) / (gamf - gamc)
375 ELSE
376 dfac = (gamma(i) - gamc) / (gamf - gamc)
377 dam(i) = dam(i) + expn * dfac**(expn-one) * dgamm(i) / (gamf - gamc)
378 END IF
379 IF (dam(i) >= one) THEN
380 dam(i) = one
381 off(i) = four_over_5
382 END IF
383 dmg(i) = dam(i)
384 END IF ! GAMMA > GAMC
385 ENDDO
386c Update damaged stress
387 DO i = 1, nel
388 dmg_scale(i) = one - dam(i)
389 END DO
390 END IF
391c-------------------------
392 soundsp(1:nel) = ssp
393c-------------------------
394 IF (nindf > 0) THEN
395 DO ii=1,nindf
396#include "lockon.inc"
397 WRITE(iout, 1000) ngl(indf(ii))
398 WRITE(istdo,1100) ngl(indf(ii)),time
399#include "lockoff.inc"
400 ENDDO
401 END IF
402c-----------------------------------------------------------------
403 1000 FORMAT(1x,'START DAMAGE IN SOLID ELEMENT NUMBER ',i10)
404 1100 FORMAT(1x,'START DAMAGE IN SOLID ELEMENT NUMBER ',i10,1x,' AT TIME :',g11.4)
405c-----------------------------------------------------------------
406 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