OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps109.F File Reference
#include "implicit_f.inc"
#include "com04_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine sigeps109 (nel, ngl, nuparam, nuvar, nvartmp, numtabl, uparam, uvar, vartmp, itable, table, jthe, time, timestep, off, rho, pla, dpla, soundsp, sigy, et, temp, epsd, dpdm, depsxx, depsyy, depszz, depsxy, depsyz, depszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, inloc, dplanl, ieos, jlag, fheat, volume, seq, lplanl, pla_nl, lepsdnl, dpdt_nl)

Function/Subroutine Documentation

◆ sigeps109()

subroutine sigeps109 ( integer nel,
integer, dimension(nel), intent(in) ngl,
integer nuparam,
integer nuvar,
integer nvartmp,
integer numtabl,
intent(in) uparam,
intent(inout) uvar,
integer, dimension(nel,nvartmp), intent(inout) vartmp,
integer, dimension(numtabl), intent(in) itable,
type(ttable), dimension(ntable) table,
integer, intent(in) jthe,
time,
timestep,
intent(in) off,
intent(in) rho,
intent(inout) pla,
intent(inout) dpla,
intent(out) soundsp,
intent(out) sigy,
intent(out) et,
intent(inout) temp,
intent(inout) epsd,
intent(in) dpdm,
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(in) dplanl,
integer, intent(in) ieos,
integer, intent(in) jlag,
intent(inout) fheat,
intent(in) volume,
intent(inout) seq,
integer, intent(in) lplanl,
dimension(nel*lplanl) pla_nl,
integer, intent(in) lepsdnl,
dimension(nel*lepsdnl) dpdt_nl )

Definition at line 34 of file sigeps109.F.

44C-----------------------------------------------
45C M o d u l e s
46C-----------------------------------------------
47 USE table_mod
48 USE interface_table_mod
49C-----------------------------------------------
50C I m p l i c i t T y p e s
51C-----------------------------------------------
52#include "implicit_f.inc"
53C-----------------------------------------------
54C C o m m o n B l o c k s
55C-----------------------------------------------
56#include "com04_c.inc"
57C-----------------------------------------------
58C D u m m y A r g u m e n t s
59C-----------------------------------------------
60 INTEGER ,INTENT(IN) :: IEOS
61 INTEGER ,INTENT(IN) :: JTHE
62 INTEGER ,INTENT(IN) :: JLAG
63 INTEGER NEL,NUPARAM,NUVAR,NVARTMP,NUMTABL,INLOC
64 INTEGER ,DIMENSION(NUMTABL),INTENT(IN) :: ITABLE
65 INTEGER ,DIMENSION(NEL) ,INTENT(IN) :: NGL
66c
67 my_real :: time,timestep
68 my_real,DIMENSION(NUPARAM) ,INTENT(IN) :: uparam
69 my_real,DIMENSION(NEL) ,INTENT(IN) :: volume
70 my_real,DIMENSION(NEL) ,INTENT(IN) :: rho,dpdm,off,
71 . depsxx,depsyy,depszz,depsxy,depsyz,depszx,
72 . sigoxx,sigoyy,sigozz,sigoxy,sigoyz,sigozx,dplanl
73 my_real ,DIMENSION(NEL) ,INTENT(OUT) :: soundsp,sigy,et,
74 . signxx,signyy,signzz,signxy,signyz,signzx
75 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: pla,dpla,epsd
76 my_real ,DIMENSION(NEL,NUVAR) ,INTENT(INOUT) :: uvar
77 INTEGER ,DIMENSION(NEL,NVARTMP) ,INTENT(INOUT) :: VARTMP
78 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: temp
79 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: fheat
80 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: seq
81 INTEGER, INTENT(IN) :: LPLANL, LEPSDNL
82 my_real, DIMENSION(NEL*LPLANL) :: pla_nl
83 my_real, DIMENSION(NEL*LEPSDNL) :: dpdt_nl
84c
85 TYPE(TTABLE), DIMENSION(NTABLE) :: TABLE
86C-----------------------------------------------
87C L o c a l V a r i a b l e s
88C-----------------------------------------------
89 INTEGER I,II,NINDX,ITER,NITER,ISMOOTH,
90 . FUNC_YLD,FUNC_TEMP,FUNC_ETA,NDIM_YLD,NDIM_TEMP,NDIM_ETA
91 INTEGER ,DIMENSION(NEL) :: INDEX
92c
93 my_real :: young,bulk,lame,g,g2,g3,nu,tref,tini,eta,ldav,cp,
94 . xrate,xscale,yscale,dtinv,j2,q2,dphi_dlam,r,dlamin,alpha,alphi
95 my_real :: asrate,ddep,dfdsig2,sig_dfdsig
96c
97 my_real, DIMENSION(NEL) ::svm,svmt,yld,yld_tref,yld_temp,
98 . sxx,syy,szz,sxy,syz,szx,sigm,
99 . fact_eta,dydx,hardp,hardr,yld_i,hardp_i,hardr_i,dxdyv,dlam,phi,
100 . ftherm,tfac,pla0,normxx,normyy,normzz,normxy,
101 . normyz,normzx,dpxx,dpyy,dpzz,dpxy,dpyz,dpzx,dpla_dlam
102 my_real, DIMENSION(NEL,3) :: xvec_eta
103 my_real, DIMENSION(NEL,4) :: xvec
104 INTEGER, DIMENSION(NEL,3) :: IPOS_ETA
105 INTEGER, DIMENSION(NEL,2) :: IPOS
106C-----------------------------------------------
107 ! VARTMP(1) latest position of PLAS in TAB_YLD function
108 ! VARTMP(2) latest position of PLAS in TAB_TEMP function
109 ! VARTMP(3) latest position of TEMP in TAB_TEMP function
110 ! VARTMP(4) latest position of TEMP in TAB_ETA function
111 ! VARTMP(5) latest position of PLAS in TAB_ETA function
112C=======================================================================
113!
114 !=========================================================================
115 !< - INITIALISATION OF COMPUTATION ON TIME STEP
116 !=========================================================================
117 !< Recovering model parameters
118 young = uparam(1) !< Young modulus
119 nu = uparam(2) !< Poisson's ratio
120 eta = uparam(3) !< Thermal work coefficient
121 tref = uparam(4) !< Reference temperature
122 tini = uparam(5) !< Initial tempareture
123 ismooth = nint(uparam(6)) !< Function interpolation flag
124 xrate = uparam(7) !< Strain rate abscissa factor for eta function
125 xscale = uparam(8) !< Strain rate abscissa factor for yld function
126 yscale = uparam(9) !< Yield function scale factor
127 g = uparam(11) !< Shear modulus
128 g2 = uparam(12) !< Shear modulus * 2
129 g3 = uparam(13) !< Shear modulus * 3
130 bulk = uparam(14) !< Bulk modulus
131 lame = uparam(15) !< Lame parameter
132 cp = uparam(20) !< Specific heat coefficient
133 asrate = uparam(21) !< Filtering pulsation for plastic strain rate
134 !< Plastic strain rate filtering parameters
135 dtinv = one / max(em20, timestep)
136 alpha = asrate*timestep
137 alphi = one-alpha
138 !< No temperature calculation inside material
139 IF (jthe /= 0) eta = zero
140!
141 !< Recovering tables and functions
142 func_yld = itable(1)
143 func_temp = itable(2)
144 func_eta = itable(3)
145 ndim_yld = table(func_yld)%NDIM
146 IF (func_temp > 0) THEN
147 ndim_temp = table(func_temp)%NDIM
148 ENDIF
149 IF (func_eta > 0) THEN
150 ndim_eta = table(func_eta)%NDIM
151 ENDIF
152!
153 !< Recovering internal variables and initializations of local variables
154 DO i = 1,nel
155 pla0(i) = pla(i) !< Initial plastic strain
156 dpla(i) = zero !< Plastic strain increment initialization
157 et(i) = one !< Hourglass stabilization variable initialization
158 hardp(i) = zero !< Hardening modulus initialization
159 ENDDO
160!
161 !< Temperature & thermal softening factor initialization
162 IF (jthe == 0) THEN
163 !< Temperature initialization
164 IF (time == zero) temp(1:nel) = tini
165 !< Thermal softening factor initialization
166 IF (eta > zero) THEN
167 !< Taylor-Quinney factor
168 IF (func_eta > 0) THEN
169 IF (inloc == 0) THEN
170 xvec_eta(1:nel,1) = epsd(1:nel) * xrate
171 ELSE
172 xvec_eta(1:nel,1) = dpdt_nl(1:nel) * xrate
173 ENDIF
174 ipos_eta(1:nel,1) = 1
175 IF (ndim_eta > 1) THEN
176 xvec_eta(1:nel,2) = temp(1:nel)
177 ipos_eta(1:nel,2) = vartmp(1:nel,4)
178 END IF
179 IF (ndim_eta > 2) THEN
180 IF (inloc == 0) THEN
181 xvec_eta(1:nel,3) = pla(1:nel)
182 ELSE
183 xvec_eta(1:nel,3) = pla_nl(1:nel)
184 ENDIF
185 ipos_eta(1:nel,3) = vartmp(1:nel,5)
186 END IF
187 CALL table_vinterp(table(func_eta),nel,nel,ipos_eta,xvec_eta,
188 . fact_eta,dxdyv)
189 IF (ndim_eta > 1) vartmp(1:nel,4) = ipos_eta(1:nel,2)
190 IF (ndim_eta > 2) vartmp(1:nel,5) = ipos_eta(1:nel,3)
191 DO i=1,nel
192 ftherm(i) = min(eta*fact_eta(i), one)
193 END DO
194 ELSE
195 ftherm(1:nel) = min(eta, one)
196 END IF
197 END IF
198 ENDIF
199!
200 !=========================================================================
201 ! - COMPUTATION OF TRIAL VALUES
202 !=========================================================================
203 DO i=1,nel
204 !< Computation of the trial stress tensor
205 ldav = (depsxx(i) + depsyy(i) + depszz(i)) * lame
206 signxx(i) = sigoxx(i) + depsxx(i)*g2 + ldav
207 signyy(i) = sigoyy(i) + depsyy(i)*g2 + ldav
208 signzz(i) = sigozz(i) + depszz(i)*g2 + ldav
209 signxy(i) = sigoxy(i) + depsxy(i)*g
210 signyz(i) = sigoyz(i) + depsyz(i)*g
211 signzx(i) = sigozx(i) + depszx(i)*g
212 !< Computation of the trace of the mean spherical stress
213 sigm(i) = (signxx(i) + signyy(i) + signzz(i)) * third
214 !< Computation of the trial deviatoric stress tensor
215 sxx(i) = signxx(i) - sigm(i)
216 syy(i) = signyy(i) - sigm(i)
217 szz(i) = signzz(i) - sigm(i)
218 sxy(i) = signxy(i)
219 syz(i) = signyz(i)
220 szx(i) = signzx(i)
221 !< Assembling Von Mises equivalent stress
222 j2 = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2) +
223 . sxy(i)**2 + syz(i)**2 + szx(i)**2
224 svm(i) = sqrt(three*j2)
225 ENDDO
226!
227 !< Computation of the initial yield stress
228 xvec(1:nel,1) = pla(1:nel)
229 xvec(1:nel,2) = epsd(1:nel) * xscale
230 ipos(1:nel,1) = vartmp(1:nel,1)
231 ipos(1:nel,2) = 1
232 CALL table2d_vinterp_log(table(func_yld),ismooth,nel,nel,ipos,
233 . xvec,yld,hardp,hardr)
234 yld(1:nel) = yld(1:nel)*yscale
235 hardp(1:nel) = hardp(1:nel)*yscale
236 vartmp(1:nel,1) = ipos(1:nel,1)
237 !< Adding temperature dependence to yield stress
238 IF (func_temp > 0) THEN
239 xvec(1:nel,2) = tref
240 ipos(1:nel,1) = vartmp(1:nel,2)
241 ipos(1:nel,2) = vartmp(1:nel,3)
242 CALL table_vinterp(table(func_temp),nel,nel,ipos,xvec,yld_tref,dydx)
243 vartmp(1:nel,2) = ipos(1:nel,1)
244 vartmp(1:nel,3) = ipos(1:nel,2)
245 xvec(1:nel,2) = temp(1:nel)
246 CALL table_vinterp(table(func_temp),nel,nel,ipos,xvec,yld_temp,dydx)
247 tfac(1:nel) = yld_temp(1:nel) / yld_tref(1:nel)
248 yld(1:nel) = yld(1:nel) * tfac(1:nel)
249 hardp(1:nel) = hardp(1:nel) * tfac(1:nel)
250 ELSE
251 tfac(1:nel) = one
252 END IF
253!
254 !=========================================================================
255 ! - COMPUTATION OF YIELD FONCTION
256 !=========================================================================
257 phi(1:nel) = svm(1:nel) - yld(1:nel)
258 nindx = 0
259 DO i=1,nel
260 IF (phi(i) >= zero .AND. off(i) == one) THEN
261 nindx = nindx + 1
262 index(nindx) = i
263 ENDIF
264 ENDDO
265!
266 !=========================================================================
267 ! - PLASTIC CORRECTION WITH CUTTING PLANE METHOD (SEMI-IMPLICIT)
268 !=========================================================================
269!
270 !< Number of iterations
271 niter = 3
272!
273 IF (nindx > 0) THEN
274!
275 !< Loop over the iterations
276 DO iter = 1,niter
277#include "vectorize.inc"
278 !< Loop over yielding elements
279 DO ii = 1, nindx
280 i = index(ii)
281!
282 ! Note: in this part, the purpose is to compute for each iteration
283 ! a plastic multiplier allowing to update internal variables to satisfy
284 ! the consistency condition using the cutting plane semi-implicit
285 ! iterative procedure.
286 ! Its expression at each iteration is : DLAMBDA = - PHI/DPHI_DLAMBDA
287 ! -> PHI : current value of yield function (known)
288 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
289 ! into account of internal variables kinetic :
290 ! plasticity, temperature and damage (to compute)
291!
292 !< 1 - Computation of the normal to the yield surface
293 !-------------------------------------------------------------------
294 normxx(i) = three_half*sxx(i)/(max(svm(i),em20))
295 normyy(i) = three_half*syy(i)/(max(svm(i),em20))
296 normzz(i) = three_half*szz(i)/(max(svm(i),em20))
297 normxy(i) = three*sxy(i)/(max(svm(i),em20))
298 normyz(i) = three*syz(i)/(max(svm(i),em20))
299 normzx(i) = three*szx(i)/(max(svm(i),em20))
300!
301 !< 2 - Computation of DPHI_DLAMBDA
302 !-------------------------------------------------------------------
303!
304 ! a) Derivative with respect stress increments tensor DSIG
305 ! ----------------------------------------------------------------
306 dfdsig2 = normxx(i) * normxx(i) * g2
307 . + normyy(i) * normyy(i) * g2
308 . + normzz(i) * normzz(i) * g2
309 . + normxy(i) * normxy(i) * g
310 . + normyz(i) * normyz(i) * g
311 . + normzx(i) * normzx(i) * g
312!
313 ! b) Derivative of dPLA with respect to DLAM
314 ! ----------------------------------------------------------------
315 sig_dfdsig = signxx(i) * normxx(i)
316 . + signyy(i) * normyy(i)
317 . + signzz(i) * normzz(i)
318 . + signxy(i) * normxy(i)
319 . + signyz(i) * normyz(i)
320 . + signzx(i) * normzx(i)
321 dpla_dlam(i) = sig_dfdsig / yld(i)
322!
323 ! c) Assemble the derivation of the yield function w.r.t. lambda
324 ! ----------------------------------------------------------------
325 dphi_dlam = - dfdsig2 + hardp(i)*dpla_dlam(i)
326 dphi_dlam = sign(max(abs(dphi_dlam),em20),dphi_dlam)
327!
328 !< 3 - Computation of the plastic multiplier
329 !-------------------------------------------------------------------
330 dlam(i) = - phi(i) / dphi_dlam
331!
332 !< 4 - Update the plastic strain related variables
333 !-------------------------------------------------------------------
334 !< Plastic strain increment on the iteration
335 ddep = dpla_dlam(i)*dlam(i)
336 !< Plastic strain increment on the time step
337 dpla(i) = max(dpla(i) + ddep,zero)
338 !< Update the plastic strain
339 pla(i) = pla0(i) + dpla(i)
340 !< Plastic strain tensor increment on the iteration
341 dpxx(i) = dlam(i)*normxx(i)
342 dpyy(i) = dlam(i)*normyy(i)
343 dpzz(i) = dlam(i)*normzz(i)
344 dpxy(i) = dlam(i)*normxy(i)
345 dpyz(i) = dlam(i)*normyz(i)
346 dpzx(i) = dlam(i)*normzx(i)
347!
348 ENDDO
349!
350 !< 5 - Update the yield stress and its derivative
351 !---------------------------------------------------------------------
352 xvec(1:nel,1:2) = zero
353 ipos(1:nel,1:2) = 0
354 DO ii = 1, nindx
355 i = index(ii)
356 xvec(ii,1) = pla(i)
357 xvec(ii,2) = epsd(i)
358 ipos(ii,1) = vartmp(i,1)
359 ipos(ii,2) = 1
360 ENDDO
361 CALL table2d_vinterp_log(table(func_yld),ismooth,nel,nindx,ipos,xvec,yld_i,hardp_i,hardr_i)
362 DO ii = 1, nindx
363 i = index(ii)
364 vartmp(i,1) = ipos(ii,1)
365 hardp(i) = hardp_i(ii)*yscale*tfac(i)
366 yld(i) = yld_i(ii)*yscale*tfac(i)
367 ENDDO
368!
369 !< 6 - Update the stress tensor and the yield function
370 !---------------------------------------------------------------------
371#include "vectorize.inc"
372 DO ii = 1, nindx
373 i = index(ii)
374!
375 !< Update the Cauchy stress tensor
376 signxx(i) = signxx(i) - dpxx(i)*g2
377 signyy(i) = signyy(i) - dpyy(i)*g2
378 signzz(i) = signzz(i) - dpzz(i)*g2
379 signxy(i) = signxy(i) - dpxy(i)*g
380 signyz(i) = signyz(i) - dpyz(i)*g
381 signzx(i) = signzx(i) - dpzx(i)*g
382!
383 !< Computation of the new trace of the mean spherical stress
384 sigm(i) = (signxx(i) + signyy(i) + signzz(i)) * third
385!
386 !< Computation of the new deviatoric stress tensor
387 sxx(i) = signxx(i) - sigm(i)
388 syy(i) = signyy(i) - sigm(i)
389 szz(i) = signzz(i) - sigm(i)
390 sxy(i) = signxy(i)
391 syz(i) = signyz(i)
392 szx(i) = signzx(i)
393!
394 !< Assembling the new Von Mises equivalent stress
395 j2 = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2) +
396 . sxy(i)**2 + syz(i)**2 + szx(i)**2
397 svm(i) = sqrt(three*j2)
398!
399 !< Update the yield function
400 phi(i) = svm(i) - yld(i)
401!
402 ENDDO
403 !< End of the loop over yielding elements
404c
405 ENDDO
406 ! End of the loop over the iterations
407 !=======================================================================
408 ! - END OF PLASTIC CORRECTION WITH CUTTING PLANE ITERATIVE METHOD
409 !=======================================================================
410!
411 !< Update the hourglass stabilization variable
412 DO ii = 1, nindx
413 i = index(ii)
414 et(i) = hardp(i) / (hardp(i) + young)
415 ENDDO
416!
417 ! Update the temperature or heat energy due to plastic work
418 IF (eta > zero .AND. inloc == 0) THEN
419 IF (jthe /= 0 .AND. jlag /= 0) THEN
420 DO ii = 1, nindx
421 i = index(ii)
422 fheat(i) = fheat(i) + ftherm(i)*yld(i)*dpla(i) * volume(i)
423 ENDDO
424 ELSE IF (cp > zero) THEN
425 DO ii = 1, nindx
426 i = index(ii)
427 temp(i) = temp(i) + ftherm(i)*yld(i)*dpla(i) / (cp*rho(i))
428 ENDDO
429 ENDIF
430 ENDIF
431!
432 ENDIF
433!
434 !< Non-local thickness variation and temperature update
435 !-------------------------------------------------------------------------
436 IF (inloc > 0 .AND. eta > zero) THEN
437 IF (jthe /= 0 .AND. jlag /= 0) THEN
438 DO i = 1, nel
439 IF (off(i) == one) THEN
440 fheat(i) = fheat(i) +
441 . ftherm(i)*yld(i)*max(dplanl(i),zero)*volume(i)
442 ENDIF
443 ENDDO
444 ELSE IF (cp > zero) THEN
445 DO i = 1,nel
446 IF (off(i) == one) THEN
447 temp(i) = temp(i) + ftherm(i)*yld(i)*max(dplanl(i),zero)/(cp*rho(i))
448 ENDIF
449 ENDDO
450 ENDIF
451 ENDIF
452!
453 !< Equation Of State (EOS) specific treatment
454 !-------------------------------------------------------------------------
455 ! if EOS is used, material law calculates only deviatoric stress tensor
456 ! sound speed depends on pressure derivative over volume change
457 ! calculated in EOS
458 IF (ieos > 0) THEN
459 signxx(1:nel) = signxx(1:nel) - sigm(1:nel)
460 signyy(1:nel) = signyy(1:nel) - sigm(1:nel)
461 signzz(1:nel) = signzz(1:nel) - sigm(1:nel)
462 signxy(1:nel) = signxy(1:nel)
463 signyz(1:nel) = signyz(1:nel)
464 signzx(1:nel) = signzx(1:nel)
465 soundsp(1:nel)= sqrt((dpdm(1:nel) + four_over_3*g) / rho(1:nel))
466 ELSE
467 soundsp(1:nel)= sqrt((bulk + four_over_3*g) / rho(1:nel))
468 END IF
469!
470 !< Storing output values
471 !-------------------------------------------------------------------------
472 DO i=1,nel
473 !< Equivalent stress
474 seq(i) = svm(i)
475 !< Yield stress
476 sigy(i) = yld(i)
477 !< Plastic strain rate
478 epsd(i) = alpha*dpla(i)*dtinv + alphi*epsd(i)
479 ENDDO
480!
my_real
#define my_real
Definition cppsort.cpp:32
alpha
#define alpha
Definition eval.h:35
interface_table_mod::table2d_vinterp_log
Definition table_mod.F:261
interface_table_mod::table_vinterp
Definition table_mod.F:251
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
interface_table_mod
Definition table_mod.F:238
table_mod
Definition table_mod.F:112