OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps109c.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!|| sigeps109c ../engine/source/materials/mat/mat109/sigeps109c.F
25!||--- called by ------------------------------------------------------
26!|| mulawc ../engine/source/materials/mat_share/mulawc.F90
27!||--- calls -----------------------------------------------------
28!|| table2d_vinterp_log ../engine/source/tools/curve/table2d_vinterp_log.F
29!|| table_vinterp ../engine/source/tools/curve/table_tools.F
30!||--- uses -----------------------------------------------------
31!|| interface_table_mod ../engine/share/modules/table_mod.F
32!|| table_mod ../engine/share/modules/table_mod.F
33!||====================================================================
34 SUBROUTINE sigeps109c(
35 1 NEL ,NGL ,NUPARAM ,NUVAR ,NVARTMP ,NUMTABL ,
36 2 UPARAM ,UVAR ,VARTMP ,ITABLE ,TABLE ,JTHE ,
37 3 TIME ,TIMESTEP,OFF ,RHO ,PLA ,DPLA ,
38 4 SOUNDSP ,SIGY ,ET ,TEMP ,EPSD ,GS ,
39 5 DEPSXX ,DEPSYY ,DEPSXY ,DEPSYZ ,DEPSZX ,
40 6 SIGOXX ,SIGOYY ,SIGOXY ,SIGOYZ ,SIGOZX ,
41 7 SIGNXX ,SIGNYY ,SIGNXY ,SIGNYZ ,SIGNZX ,
42 8 THK ,THKLY ,INLOC ,DPLANL ,LOFF ,SEQ ,
43 9 LPLANL ,PLA_NL ,LEPSDNL ,DPDT_NL )
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 NEL,NUPARAM,NUVAR,NVARTMP,NUMTABL,JTHE,INLOC
61 INTEGER ,DIMENSION(NUMTABL),INTENT(IN) :: ITABLE
62 INTEGER ,DIMENSION(NEL) ,INTENT(IN) :: NGL
63c
64 my_real :: TIME,TIMESTEP
65 my_real,DIMENSION(NUPARAM) ,INTENT(IN) :: UPARAM
66 my_real,DIMENSION(NEL) ,INTENT(IN) :: RHO,OFF,GS,THKLY,
67 . DEPSXX,DEPSYY,DEPSXY,DEPSYZ,DEPSZX,
68 . SIGOXX,SIGOYY,SIGOXY,SIGOYZ,SIGOZX,DPLANL
69 my_real ,DIMENSION(NEL) ,INTENT(OUT) :: soundsp,sigy,et,
70 . signxx,signyy,signxy,signyz,signzx
71 my_real ,DIMENSION(NEL) ,INTENT(INOUT) :: pla,dpla,epsd,temp,thk,seq
72 my_real ,DIMENSION(NEL,NUVAR) ,INTENT(INOUT) :: uvar
73 INTEGER ,DIMENSION(NEL,NVARTMP) ,INTENT(INOUT) :: VARTMP
74c
75 TYPE(ttable), DIMENSION(NTABLE) :: TABLE
76 my_real, DIMENSION(NEL), INTENT(IN) :: LOFF
77 INTEGER, INTENT(IN) :: LPLANL, LEPSDNL
78 my_real, DIMENSION(NEL*LPLANL) , INTENT(IN) :: pla_nl
79 my_real, DIMENSION(NEL*LEPSDNL), INTENT(IN) :: dpdt_nl
80C-----------------------------------------------
81C L o c a l V a r i a b l e s
82C-----------------------------------------------
83 INTEGER I,II,NINDX,ITER,NITER,ISMOOTH,
84 . func_yld,func_temp,func_eta,ndim_yld,ndim_temp,ndim_eta
85 INTEGER ,DIMENSION(NEL) :: INDEX
86c
87 my_real :: young,g,a11,a12,nu,tref,tini,eta,
88 . xrate,xscale,yscale,dtinv,j2,q2,dphi_dlam,alpha,alphi,
89 . r,dfdsig2,sig_dfdsig,ddep,asrate
90c
91 my_real, DIMENSION(NEL) ::svm,svmt,yld,yld_tref,yld_temp,
92 . sxx,syy,sxy,syz,szx,sigm,stxx,styy,stxy,styz,stzx,
93 . fact_eta,dydx,hardp,hardr,yld_i,hardp_i,hardr_i,dxdyv,dlam,phi,
94 . ftherm,tfac,depszz,normxx,normyy,normxy,depspxx,depspyy,depspxy,
95 . stzz,szz,dpxx,dpyy,dpxy,dpla_dlam,pla0,dezz
96 my_real, DIMENSION(NEL,3) :: xvec_eta
97 my_real, DIMENSION(NEL,4) :: xvec
98 INTEGER, DIMENSION(NEL,3) :: IPOS_ETA
99 INTEGER, DIMENSION(NEL,2) :: IPOS
100!-------------------------------------------------------------------------------
101 ! VARTMP(1) latest position of PLAS in TAB_YLD function
102 ! VARTMP(2) latest position of PLAS in TAB_TEMP function
103 ! VARTMP(3) latest position of TEMP in TAB_TEMP function
104 ! VARTMP(4) latest position of TEMP in TAB_ETA function
105 ! VARTMP(5) latest position of PLAS in TAB_ETA function
106!===============================================================================
107!
108 !=========================================================================
109 !< - INITIALISATION OF COMPUTATION ON TIME STEP
110 !=========================================================================
111 !< Recovering model parameters
112 young = uparam(1) !< Young modulus
113 nu = uparam(2) !< Poisson's ratio
114 eta = uparam(3) !< Thermal work coefficient
115 tref = uparam(4) !< Reference temperature
116 tini = uparam(5) !< Initial tempareture
117 ismooth = nint(uparam(6)) !< Function interpolation flag
118 xrate = uparam(7) !< Strain rate abscissa factor for eta function
119 xscale = uparam(8) !< Strain rate abscissa factor for yld function
120 yscale = uparam(9) !< Yield function scale factor
121 g = uparam(11) !< Shear modulus
122 a11 = uparam(16) !< 2D Elastic matrix diagonal coefficient
123 a12 = uparam(17) !< 2D Elastic matrix off-diagonal coefficient
124 asrate = uparam(21) !< Filtering pulsation for plastic strain rate
125 !< Plastic strain rate filtering parameters
126 dtinv = one / max(em20, timestep)
127 alpha = asrate*timestep
128 alphi = one-alpha
129 !< No temperature calculation inside material
130 IF (jthe /= 0) eta = zero
131!
132 !< Recovering tables and functions
133 func_yld = itable(1)
134 func_temp = itable(2)
135 func_eta = itable(3)
136 ndim_yld = table(func_yld)%NDIM
137 IF (func_temp > 0) THEN
138 ndim_temp = table(func_temp)%NDIM
139 ENDIF
140 IF (func_eta > 0) THEN
141 ndim_eta = table(func_eta)%NDIM
142 ENDIF
143!
144 !< Recovering internal variables and initializations of local variables
145 DO i = 1,nel
146 pla0(i) = pla(i) !< Initial plastic strain
147 dpla(i) = zero !< Plastic strain increment initialization
148 et(i) = one !< Hourglass stabilization variable initialization
149 hardp(i) = zero !< Hardening modulus initialization
150 dezz(i) = zero !< Thickness variation initialization
151 ENDDO
152!
153 !< Temperature & thermal softening factor initialization
154 IF (jthe == 0) THEN
155 !< Temperature initialization
156 IF (time == zero) temp(1:nel) = tini
157 !< Thermal softening factor initialization
158 IF (eta > zero) THEN
159 !< Taylor-Quinney factor
160 IF (func_eta > 0) THEN
161 IF (inloc == 0) THEN
162 xvec_eta(1:nel,1) = epsd(1:nel) * xrate
163 ELSE
164 xvec_eta(1:nel,1) = dpdt_nl(1:nel) * xrate
165 ENDIF
166 ipos_eta(1:nel,1) = 1
167 IF (ndim_eta > 1) THEN
168 xvec_eta(1:nel,2) = temp(1:nel)
169 ipos_eta(1:nel,2) = vartmp(1:nel,4)
170 END IF
171 IF (ndim_eta > 2) THEN
172 IF (inloc == 0) THEN
173 xvec_eta(1:nel,3) = pla(1:nel)
174 ELSE
175 xvec_eta(1:nel,3) = pla_nl(1:nel)
176 ENDIF
177 ipos_eta(1:nel,3) = vartmp(1:nel,5)
178 END IF
179 CALL table_vinterp(table(func_eta),nel,nel,ipos_eta,xvec_eta,
180 . fact_eta,dxdyv)
181 IF (ndim_eta > 1) vartmp(1:nel,4) = ipos_eta(1:nel,2)
182 IF (ndim_eta > 2) vartmp(1:nel,5) = ipos_eta(1:nel,3)
183 DO i=1,nel
184 ftherm(i) = min(eta*fact_eta(i), one)
185 END DO
186 ELSE
187 ftherm(1:nel) = min(eta, one)
188 END IF
189 END IF
190 ENDIF
191!
192 !=========================================================================
193 ! - COMPUTATION OF TRIAL VALUES
194 !=========================================================================
195 DO i=1,nel
196 !< Computation of the trial stress tensor
197 signxx(i) = sigoxx(i) + a11*depsxx(i) + a12*depsyy(i)
198 signyy(i) = sigoyy(i) + a11*depsyy(i) + a12*depsxx(i)
199 signxy(i) = sigoxy(i) + depsxy(i)*g
200 signyz(i) = sigoyz(i) + depsyz(i)*gs(i)
201 signzx(i) = sigozx(i) + depszx(i)*gs(i)
202 !< Computation of the trace of the mean spherical stress
203 sigm(i) = (signxx(i) + signyy(i)) * third
204 !< Computation of the trial deviatoric stress tensor
205 sxx(i) = signxx(i) - sigm(i)
206 syy(i) = signyy(i) - sigm(i)
207 szz(i) = -sigm(i)
208 sxy(i) = signxy(i)
209 !< Assembling Von Mises equivalent stress
210 j2 = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 + two*sxy(i)**2)
211 svm(i) = sqrt(three*j2)
212 ENDDO
213!
214 !< Computation of the initial yield stress
215 xvec(1:nel,1) = pla(1:nel)
216 xvec(1:nel,2) = epsd(1:nel) * xscale
217 ipos(1:nel,1) = vartmp(1:nel,1)
218 ipos(1:nel,2) = 1
219 CALL table2d_vinterp_log(table(func_yld),ismooth,nel,nel,ipos,
220 . xvec,yld,hardp,hardr)
221 yld(1:nel) = yld(1:nel)*yscale
222 hardp(1:nel) = hardp(1:nel)*yscale
223 vartmp(1:nel,1) = ipos(1:nel,1)
224 !< Adding temperature dependence to yield stress
225 IF (func_temp > 0) THEN
226 xvec(1:nel,2) = tref
227 ipos(1:nel,1) = vartmp(1:nel,2)
228 ipos(1:nel,2) = vartmp(1:nel,3)
229 CALL table_vinterp(table(func_temp),nel,nel,ipos,xvec,yld_tref,dydx)
230 vartmp(1:nel,2) = ipos(1:nel,1)
231 vartmp(1:nel,3) = ipos(1:nel,2)
232 xvec(1:nel,2) = temp(1:nel)
233 CALL table_vinterp(table(func_temp),nel,nel,ipos,xvec,yld_temp,dydx)
234 tfac(1:nel) = yld_temp(1:nel) / yld_tref(1:nel)
235 yld(1:nel) = yld(1:nel) * tfac(1:nel)
236 hardp(1:nel) = hardp(1:nel) * tfac(1:nel)
237 ELSE
238 tfac(1:nel) = one
239 END IF
240!
241 !=========================================================================
242 ! - COMPUTATION OF YIELD FONCTION
243 !=========================================================================
244 phi(1:nel) = svm(1:nel) - yld(1:nel)
245 nindx = 0
246 DO i=1,nel
247 IF (phi(i) >= zero .AND. off(i) == one) THEN
248 nindx = nindx + 1
249 index(nindx) = i
250 ENDIF
251 ENDDO
252!
253 !=========================================================================
254 ! - PLASTIC CORRECTION WITH CUTTING PLANE METHOD (SEMI-IMPLICIT)
255 !=========================================================================
256!
257 !< Number of iterations
258 niter = 3
259!
260 IF (nindx > 0) THEN
261!
262 !< Loop over the iterations
263 DO iter = 1,niter
264#include "vectorize.inc"
265 !< Loop over yielding elements
266 DO ii = 1, nindx
267 i = index(ii)
268!
269 ! Note: in this part, the purpose is to compute for each iteration
270 ! a plastic multiplier allowing to update internal variables to satisfy
271 ! the consistency condition using the cutting plane semi-implicit
272 ! iterative procedure.
273 ! Its expression at each iteration is : DLAMBDA = - PHI/DPHI_DLAMBDA
274 ! -> PHI : current value of yield function (known)
275 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
276 ! into account of internal variables kinetic :
277 ! plasticity, temperature and damage (to compute)
278!
279 !< 1 - Computation of the normal to the yield surface
280 !-------------------------------------------------------------------
281 normxx(i) = three_half*sxx(i)/(max(svm(i),em20))
282 normyy(i) = three_half*syy(i)/(max(svm(i),em20))
283 normxy(i) = three*sxy(i)/(max(svm(i),em20))
284!
285 !< 2 - Computation of DPHI_DLAMBDA
286 !-------------------------------------------------------------------
287!
288 ! a) Derivative with respect stress increments tensor DSIG
289 ! ----------------------------------------------------------------
290 dfdsig2 = normxx(i) * (a11*normxx(i) + a12*normyy(i))
291 . + normyy(i) * (a11*normyy(i) + a12*normxx(i))
292 . + normxy(i) * normxy(i) * g
293!
294 ! b) Derivative of dPLA with respect to DLAM
295 ! ----------------------------------------------------------------
296 sig_dfdsig = signxx(i) * normxx(i)
297 . + signyy(i) * normyy(i)
298 . + signxy(i) * normxy(i)
299 dpla_dlam(i) = sig_dfdsig / yld(i)
300!
301 ! c) Assemble the derivation of the yield function w.r.t. lambda
302 ! ----------------------------------------------------------------
303 dphi_dlam = - dfdsig2 + hardp(i)*dpla_dlam(i)
304 dphi_dlam = sign(max(abs(dphi_dlam),em20),dphi_dlam)
305!
306 !< 3 - Computation of the plastic multiplier
307 !-------------------------------------------------------------------
308 dlam(i) = - phi(i) / dphi_dlam
309!
310 !< 4 - Update the plastic strain related variables
311 !-------------------------------------------------------------------
312 !< Plastic strain increment on the iteration
313 ddep = dpla_dlam(i)*dlam(i)
314 !< Plastic strain increment on the time step
315 dpla(i) = max(dpla(i) + ddep,zero)
316 !< Update the plastic strain
317 pla(i) = pla0(i) + dpla(i)
318 !< Plastic strain tensor increment on the iteration
319 dpxx(i) = dlam(i)*normxx(i)
320 dpyy(i) = dlam(i)*normyy(i)
321 dpxy(i) = dlam(i)*normxy(i)
322!
323 ENDDO
324!
325 !< 5 - Update the yield stress and its derivative
326 !---------------------------------------------------------------------
327 xvec(1:nel,1:2) = zero
328 ipos(1:nel,1:2) = 0
329 DO ii = 1, nindx
330 i = index(ii)
331 xvec(ii,1) = pla(i)
332 xvec(ii,2) = epsd(i)
333 ipos(ii,1) = vartmp(i,1)
334 ipos(ii,2) = 1
335 ENDDO
336 CALL table2d_vinterp_log(table(func_yld),ismooth,nel,nindx,ipos,xvec,yld_i,hardp_i,hardr_i)
337 DO ii = 1, nindx
338 i = index(ii)
339 vartmp(i,1) = ipos(ii,1)
340 hardp(i) = hardp_i(ii)*yscale*tfac(i)
341 yld(i) = yld_i(ii)*yscale*tfac(i)
342 ENDDO
343!
344 !< 6 - Update the stress tensor and the yield function
345 !---------------------------------------------------------------------
346#include "vectorize.inc"
347 DO ii = 1, nindx
348 i = index(ii)
349!
350 !< Update the Cauchy stress tensor
351 signxx(i) = signxx(i) - a11*dpxx(i) - a12*dpyy(i)
352 signyy(i) = signyy(i) - a12*dpxx(i) - a11*dpyy(i)
353 signxy(i) = signxy(i) - g*dpxy(i)
354!
355 !< Computation of the new trace of the mean spherical stress
356 sigm(i) = (signxx(i) + signyy(i)) * third
357!
358 !< Computation of the new deviatoric stress tensor
359 sxx(i) = signxx(i) - sigm(i)
360 syy(i) = signyy(i) - sigm(i)
361 szz(i) = -sigm(i)
362 sxy(i) = signxy(i)
363!
364 !< Assembling the new Von Mises equivalent stress
365 j2 = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 + two*sxy(i)**2)
366 svm(i) = sqrt(three*j2)
367!
368 !< Update the yield function
369 phi(i) = svm(i) - yld(i)
370!
371 !< Update the thickness variation
372 IF (inloc == 0) dezz(i) = dezz(i) - dpxx(i) - dpyy(i)
373!
374 ENDDO
375 !< End of the loop over yielding elements
376 ENDDO
377 ! End of the loop over the iterations
378 !=======================================================================
379 ! - END OF PLASTIC CORRECTION WITH CUTTING PLANE ITERATIVE METHOD
380 !=======================================================================
381!
382 !< Update the hourglass stabilization variable
383 DO ii = 1, nindx
384 i = index(ii)
385 et(i) = hardp(i) / (hardp(i) + young)
386 ENDDO
387!
388 !< Update the temperature
389 IF (jthe == 0 .AND. eta > zero .AND. inloc == 0) THEN
390 DO ii = 1, nindx
391 i = index(ii)
392 temp(i) = temp(i) + ftherm(i)*yld(i)*dpla(i)
393 ENDDO
394 ENDIF
395 ENDIF
396!
397 !< Non-local thickness variation and temperature update
398 !-------------------------------------------------------------------------
399 IF (inloc > 0) THEN
400 DO i = 1,nel
401 IF (loff(i) == one) THEN
402 dezz(i) = max(dplanl(i),zero)*half*(signxx(i)+signyy(i))/max(yld(i),em20)
403 dezz(i) = - nu*((signxx(i)-sigoxx(i)+signyy(i)-sigoyy(i))/young) - dezz(i)
404 IF (jthe == 0 .AND. eta > zero) THEN
405 temp(i) = temp(i) + ftherm(i)*yld(i)*max(dplanl(i),zero)
406 ENDIF
407 ENDIF
408 ENDDO
409 ENDIF
410!
411 !< Storing output values and update soundspeed and thickness
412 !-------------------------------------------------------------------------
413 DO i=1,nel
414 !< Equivalent stress
415 seq(i) = svm(i)
416 !< Sound speed
417 soundsp(i) = sqrt(a11/rho(i))
418 !< Yield stress
419 sigy(i) = yld(i)
420 !< Plastic strain rate
421 epsd(i) = alpha*dpla(i)*dtinv + alphi*epsd(i)
422 !< Thickness variation
423 IF (inloc > 0) THEN
424 IF (loff(i) == one) THEN
425 dezz(i) = -nu*(signxx(i)-sigoxx(i)+signyy(i)-sigoyy(i))/young
426 dezz(i) = dezz(i) - max(dplanl(i),zero)*half*(signxx(i)+signyy(i))/max(yld(i),em20)
427 ENDIF
428 ELSE
429 dezz(i) = -nu*(signxx(i)-sigoxx(i)+signyy(i)-sigoyy(i))/young + dezz(i)
430 ENDIF
431 thk(i) = thk(i) + dezz(i)*thkly(i)*off(i)
432 ENDDO
433!
434 END SUBROUTINE sigeps109c
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
sigeps109c
subroutine sigeps109c(nel, ngl, nuparam, nuvar, nvartmp, numtabl, uparam, uvar, vartmp, itable, table, jthe, time, timestep, off, rho, pla, dpla, soundsp, sigy, et, temp, epsd, gs, depsxx, depsyy, depsxy, depsyz, depszx, sigoxx, sigoyy, sigoxy, sigoyz, sigozx, signxx, signyy, signxy, signyz, signzx, thk, thkly, inloc, dplanl, loff, seq, lplanl, pla_nl, lepsdnl, dpdt_nl)
Definition sigeps109c.F:44
table_mod::ttable
Definition table_mod.F:125