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