OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps190.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!|| sigeps190 ../engine/source/materials/mat/mat190/sigeps190.F
25!||--- called by ------------------------------------------------------
26!|| mulaw ../engine/source/materials/mat_share/mulaw.F90
27!||--- calls -----------------------------------------------------
28!|| condamage ../engine/source/materials/mat/mat190/condamage.F
29!|| conversion ../engine/source/materials/mat/mat190/conversion.F
30!|| prodata ../engine/source/materials/tools/prodATA.F
31!|| valpvec_v ../engine/source/materials/mat/mat033/sigeps33.F
32!|| valpvecdp_v ../engine/source/materials/mat/mat033/sigeps33.F
33!||--- uses -----------------------------------------------------
34!|| condamage_mod ../engine/source/materials/mat/mat190/condamage.F
35!|| conversion_mod ../engine/source/materials/mat/mat190/conversion.F
36!|| matparam_def_mod ../common_source/modules/mat_elem/matparam_def_mod.F90
37!||====================================================================
38 SUBROUTINE sigeps190(
39 1 NEL ,NUVAR ,RHO ,ET ,
40 2 EPSPXX ,EPSPYY ,EPSPZZ ,EPSPXY ,EPSPYZ ,EPSPZX ,
41 3 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
42 4 MFXX ,MFXY ,MFXZ ,MFYX ,MFYY ,MFYZ ,
43 5 MFZX ,MFZY ,MFZZ ,SOUNDSP ,VISCMAX ,UVAR ,
44 7 NUMTABL,MATPARAM,NVARTMP ,VARTMP )
45C-----------------------------------------------
46C M o d u l e s
47C-----------------------------------------------
48 USE conversion_mod
49 USE condamage_mod
50 USE matparam_def_mod
51C-----------------------------------------------
52C Material law for isotropic path dependent recoverable foam
53C-----------------------------------------------
54C I m p l i c i t T y p e s
55C-----------------------------------------------
56#include "implicit_f.inc"
57C-----------------------------------------------
58#include "comlock.inc"
59#include "scr05_c.inc"
60#include "impl1_c.inc"
61#include "com01_c.inc"
62#include "mvsiz_p.inc"
63C-----------------------------------------------
64C G l o b a l P a r a m e t e r s
65C-----------------------------------------------
66C---------+---------+---+---+--------------------------------------------
67C VAR | SIZE |TYP| RW| DEFINITION
68C---------+---------+---+---+--------------------------------------------
69C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
70C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
71C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
72C---------+---------+---+---+--------------------------------------------
73C---------+---------+---+---+--------------------------------------------
74C TIME | 1 | F | R | CURRENT TIME
75C TIMESTEP| 1 | F | R | CURRENT TIME STEP
76C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
77C RHO0 | NEL | F | R | INITIAL DENSITY
78C RHO | NEL | F | R | DENSITY
79C VOLUME | NEL | F | R | VOLUME
80C EINT | NEL | F | R | TOTAL INTERNAL ENERGY
81C EPSPXX | NEL | F | R | STRAIN RATE XX
82C EPSPYY | NEL | F | R | STRAIN RATE YY
83C ... | | | |
84C---------+---------+---+---+--------------------------------------------
85C SIGNXX | NEL | F | W | NEW ELASTO PLASTIC STRESS XX
86C SIGNYY | NEL | F | W | NEW ELASTO PLASTIC STRESS YY
87C ... | | | |
88C SIGVXX | NEL | F | W | VISCOUS STRESS XX
89C SIGVYY | NEL | F | W | VISCOUS STRESS YY
90C ... | | | |
91C SOUNDSP | NEL | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
92C VISCMAX | NEL | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
93C---------+---------+---+---+--------------------------------------------
94C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
95C OFF | NEL | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
96C---------+---------+---+---+--------------------------------------------
97C-------------------------------------------------------------------------
98C-----------------------------------------------
99C I N P U T A r g u m e n t s
100C-----------------------------------------------
101 INTEGER ,INTENT(IN) :: NEL,NUVAR,NVARTMP,NUMTABL
102 !INTEGER, DIMENSION(NEL),INTENT(IN) :: NGL
103 my_real ,DIMENSION(NEL) ,INTENT(IN) :: RHO,
104 . EPSPXX,EPSPYY,EPSPZZ,EPSPXY,EPSPYZ,EPSPZX,
105 . MFXX ,MFXY ,MFXZ ,MFYX, MFYY , MFYZ ,
106 . MFZX ,MFZY ,MFZZ
107C-----------------------------------------------
108C O U T P U T A r g u m e n t s
109C-----------------------------------------------
110 my_real ,DIMENSION(NEL) ,INTENT(OUT) ::
111 . signxx ,signyy ,signzz ,signxy ,signyz ,signzx,viscmax
112C-----------------------------------------------
113C I N P U T O U T P U T A r g u m e n t s
114C-----------------------------------------------
115 my_real, DIMENSION(NEL) ,INTENT(INOUT) :: soundsp, et
116 my_real ,DIMENSION(NEL,NUVAR) ,INTENT(INOUT) :: uvar
117 INTEGER ,DIMENSION(NEL,NVARTMP) ,INTENT(INOUT) :: VARTMP
118 TYPE(MATPARAM_STRUCT_) :: MATPARAM
119C----------------------------------------------------------------
120C L O C A L V A R I B L E S
121C----------------------------------------------------------------
122 INTEGER I,IK,IJ
123 my_real, DIMENSION(NEL) ::
124 . jac ,jeq ,b11 ,b22 ,b33 ,b12 ,
125 . b23 ,b13 ,bsqr11 ,bsqr22 ,i1bar ,i2bar ,
126 . bsqr33 ,bsqr12 ,bsqr23 ,bsqr13 ,h1 ,h2 ,h3 ,
127 . zxx ,zxy ,zyz ,zxz , zyx , zzy ,
128 . zzx ,zyy ,zzz ,whysmax ,
129 . zxxold ,zxyold ,zyzold ,zxzold ,
130 . zyxold , zzyold ,zzxold ,zyyold ,zzzold ,
131 . zxxmid ,zxymid ,zyzmid ,zxzmid ,
132 . zyxmid , zzymid ,zzxmid ,zyymid ,zzzmid ,
133 . dzxx ,dzxy ,dzyz ,dzxz ,
134 . dzyx , dzzy ,dzzx ,dzyy ,dzzz ,
135 . s1 ,s2 ,s3 ,s4 ,s5 ,s6 ,
136 . fs11 ,fs12 ,fs13 ,fs21 ,fs22 ,fs23 ,
137 . fs31 ,fs32 ,fs33 ,slopemax , damage ,
138 . zxxd ,zyyd ,zzzd ,zxyd ,zyzd ,zzxd ,epst ,yld
139c
140 my_real
141 . g,k,nu,hu, shape ,gs,ks,stiffmax,aa, de,e_new,epss, e_old,
142 . ep2,ep3,ep4,ep5,ep6,ert11,ert12,ert13,ert21,e,emod,ep1,
143 . ert22,ert23,ert31,ert32,ert33,scal,
144 . epsp(3),xscal,
145C
146 . strain(nel,3),slope(nel,3),drate(nel,3),ev(nel,3),
147 . f(nel,3,3),fo(nel,3,3),c(nel,3,3), egl(nel,6),
148C
149 . spknodam(nel,6),spknorate(nel,6),spk(nel,6),sig(nel,6) ,
150 . dfirst(nel,3,6),cijkl(nel,6,6),
151c
152 . evv(mvsiz,3), val(mvsiz,3),av(mvsiz,6),dirprv(mvsiz,3,3)
153c
154c=========================================================================
155c Isotropic path dependent recoverable foam law
156c-------------------------------------------------------------------------
157!
158 !=======================================================================
159 ! - INITIALISATION OF COMPUTATION ON TIME STEP
160 !=======================================================================
161 ! Recovering model parameters
162 emod = matparam%UPARAM(1)
163 k = matparam%UPARAM(2)
164 g = matparam%UPARAM(3)
165 xscal = matparam%UPARAM(4)
166 scal = matparam%UPARAM(5)
167 hu = matparam%UPARAM(6)
168 shape = matparam%UPARAM(7)
169 stiffmax = emod*hundred
170 ! Recovering user variables
171 whysmax(1:nel) = uvar(1:nel,13)
172!
173 !=======================================================================
174 ! - COMPUTATION OF THE DEFORMATION GRADIENT
175 !=======================================================================
176 DO i=1,nel
177 f(i,1,1) = one+mfxx(i)
178 f(i,2,2) = one+mfyy(i)
179 f(i,3,3) = one+mfzz(i)
180 f(i,1,2) = mfxy(i)
181 f(i,2,3) = mfyz(i)
182 f(i,1,3) = mfxz(i)
183 f(i,2,1) = mfyx(i)
184 f(i,3,2) = mfzy(i)
185 f(i,3,1) = mfzx(i)
186 ENDDO
187!
188 !=======================================================================
189 ! - COMPUTATION OF THE GREEN-LAGRANGE TENSORS (FT * F - I)/2
190 !=======================================================================
191 CALL prodata(f, c, egl, nel)
192 DO i=1,nel
193 ! Current Green-Lagrange tensor
194 zxx(i) = egl(i,1)
195 zyy(i) = egl(i,2)
196 zzz(i) = egl(i,3)
197 zxy(i) = egl(i,4)
198 zyz(i) = egl(i,5)
199 zzx(i) = egl(i,6)
200C
201 ! Recover old Green-Lagrange tensor
202 zxxold(i) = uvar(i,1)
203 zyyold(i) = uvar(i,2)
204 zzzold(i) = uvar(i,3)
205 zxyold(i) = uvar(i,4)
206 zyzold(i) = uvar(i,5)
207 zzxold(i) = uvar(i,6)
208c
209 ! Save current Green-Lagrange tensor
210 uvar(i,1) = zxx(i)
211 uvar(i,2) = zyy(i)
212 uvar(i,3) = zzz(i)
213 uvar(i,4) = zxy(i)
214 uvar(i,5) = zyz(i)
215 uvar(i,6) = zzx(i)
216c
217 ! Midpoint Green-Lagrange tensor at t = n + 1/2
218 zxxmid(i) = (zxx(i) + zxxold(i))/two
219 zyymid(i) = (zyy(i) + zyyold(i))/two
220 zzzmid(i) = (zzz(i) + zzzold(i))/two
221 zxymid(i) = (zxy(i) + zxyold(i))/two
222 zyzmid(i) = (zyz(i) + zyzold(i))/two
223 zzxmid(i) = (zzx(i) + zzxold(i))/two
224c
225 ! Deformation Jacobian
226 jac(i) = f(i,1,1)*f(i,2,2)*f(i,3,3) - f(i,1,1)*f(i,2,3)*f(i,3,2) -
227 . f(i,3,3)*f(i,1,2)*f(i,2,1) + f(i,1,2)*f(i,2,3)*f(i,3,1) +
228 . f(i,2,1)*f(i,3,2)*f(i,1,3) - f(i,2,2)*f(i,3,1)*f(i,1,3)
229 jac(i) = max(jac(i),one)
230c
231 ENDDO
232c
233 !=======================================================================
234 ! - COMPUTATION OF THE GREEN-LAGRANGE EIGENVALUES AND VECTORS
235 !=======================================================================
236 DO i=1,nel
237 av(i,1) = zxxmid(i)
238 av(i,2) = zyymid(i)
239 av(i,3) = zzzmid(i)
240 av(i,4) = zxymid(i)
241 av(i,5) = zyzmid(i)
242 av(i,6) = zzxmid(i)
243 ENDDO
244 IF (iresp == 1) THEN
245 CALL valpvecdp_v(av,evv,dirprv,nel)
246 ELSE
247 CALL valpvec_v(av,evv,dirprv,nel)
248 ENDIF
249c
250 !=======================================================================
251 ! - COMPUTATION OF THE STRAIN RATES IN PRINCIPAL DIRECTIONS
252 !=======================================================================
253 DO i=1,nel
254c
255 ! True strain rates
256 ep1 = epspxx(i)
257 ep2 = epspyy(i)
258 ep3 = epspzz(i)
259 ep4 = half*epspxy(i)
260 ep5 = half*epspyz(i)
261 ep6 = half*epspzx(i)
262c
263 ! Rotation in principal directions
264 ert11 = dirprv(i,1,1)*ep1 + dirprv(i,2,1)*ep4 + dirprv(i,3,1)*ep6
265 ert12 = dirprv(i,1,2)*ep1 + dirprv(i,2,2)*ep4 + dirprv(i,3,2)*ep6
266 ert13 = dirprv(i,1,3)*ep1 + dirprv(i,2,3)*ep4 + dirprv(i,3,3)*ep6
267 ert21 = dirprv(i,1,1)*ep4 + dirprv(i,2,1)*ep2 + dirprv(i,3,1)*ep5
268 ert22 = dirprv(i,1,2)*ep4 + dirprv(i,2,2)*ep2 + dirprv(i,3,2)*ep5
269 ert23 = dirprv(i,1,3)*ep4 + dirprv(i,2,3)*ep2 + dirprv(i,3,3)*ep5
270 ert31 = dirprv(i,1,1)*ep6 + dirprv(i,2,1)*ep5 + dirprv(i,3,1)*ep3
271 ert32 = dirprv(i,1,2)*ep6 + dirprv(i,2,2)*ep5 + dirprv(i,3,2)*ep3
272 ert33 = dirprv(i,1,3)*ep6 + dirprv(i,2,3)*ep5 + dirprv(i,3,3)*ep3
273 epsp(1) = dirprv(i,1,1)*ert11 + dirprv(i,2,1)*ert21
274 . + dirprv(i,3,1)*ert31
275 epsp(2) = dirprv(i,1,2)*ert12 + dirprv(i,2,2)*ert22
276 . + dirprv(i,3,2)*ert32
277 epsp(3) = dirprv(i,1,3)*ert13 + dirprv(i,2,3)*ert23
278 . + dirprv(i,3,3)*ert33
279c
280 !(e > 0 compression and e < 0 traction)
281 ev(i,1) = sqrt(two*evv(i,1) + one)
282 ev(i,2) = sqrt(two*evv(i,2) + one)
283 ev(i,3) = sqrt(two*evv(i,3) + one)
284 strain(i,1) = one - ev(i,1)
285 strain(i,2) = one - ev(i,2)
286 strain(i,3) = one - ev(i,3)
287 epst(i) = sqrt(strain(i,1)**2 + strain(i,2)**2 + strain(i,3)**2)
288 drate(i,1) = epsp(1)*(one - strain(i,1)) ! eng
289 drate(i,2) = epsp(2)*(one - strain(i,2))
290 drate(i,3) = epsp(3)*(one - strain(i,3))
291 ENDDO
292c
293 !=======================================================================
294 ! - COMPUTATION OF THE GLOBAL INCREMENTAL STIFFNESS MATRIX
295 ! (Vérifier si effet de Poisson purement numérique)
296 !=======================================================================
297 cijkl = zero
298 CALL conversion(zxxmid ,zyymid ,zzzmid ,zxymid ,zyzmid ,zzxmid ,
299 . cijkl ,dfirst ,drate ,xscal ,scal ,nel ,
300 . jac ,slope ,numtabl ,matparam%TABLE ,nvartmp ,
301 . vartmp ,yld ,matparam%TABLE(1)%NDIM )
302C
303 DO i = 1,nel
304C
305C increment of Green Lagrange strains
306C
307 dzxx(i) = zxx(i)-zxxold(i)
308 dzyy(i) = zyy(i)-zyyold(i)
309 dzzz(i) = zzz(i)-zzzold(i)
310 dzxy(i) = zxy(i)-zxyold(i)
311 dzyz(i) = zyz(i)-zyzold(i)
312 dzzx(i) = zzx(i)-zzxold(i)
313C
314C USE GAMMAS INSTEAD OF EPSILONS FOR INCREMENTS
315C
316 dzxy(i) = dzxy(i)*two
317 dzyz(i) = dzyz(i)*two
318 dzzx(i) = dzzx(i)*two
319C
320C INCREMENT HYPERELASTIC 2PK STRESSES DS=C*DE
321C
322 DO ik=1,6
323 spknorate(i,ik)= uvar(i,6+ik)
324 . + cijkl(i,1,ik)*dzxx(i)+cijkl(i,2,ik)*dzyy(i)
325 . + cijkl(i,3,ik)*dzzz(i)+cijkl(i,4,ik)*dzxy(i)
326 . + cijkl(i,5,ik)*dzyz(i)+cijkl(i,6,ik)*dzzx(i)
327 ENDDO
328 uvar(i, 7) = spknorate(i,1)
329 uvar(i, 8) = spknorate(i,2)
330 uvar(i, 9) = spknorate(i,3)
331 uvar(i,10) = spknorate(i,4)
332 uvar(i,11) = spknorate(i,5)
333 uvar(i,12) = spknorate(i,6)
334
335c add the viscous 2PK stresses ( dynamic overstress )
336c drate = dynamic principal 2PK - static principal 2PK
337C
338
339 DO ik=1,6
340 IF (ik < 4) THEN
341 spknodam(i,ik)=spknorate(i,ik)
342 . +dfirst(i,1,ik)*drate(i,1)
343 . +dfirst(i,2,ik)*drate(i,2)
344 . +dfirst(i,3,ik)*drate(i,3)
345 ELSE
346 spknodam(i,ik)=spknorate(i,ik)
347 . + (dfirst(i,1,ik)*drate(i,1)
348 . + dfirst(i,2,ik)*drate(i,2)
349 . + dfirst(i,3,ik)*drate(i,3) )*half
350 ENDIF
351 ENDDO
352c
353 ENDDO !nel
354C
355C compute damage at t(n+1)
356 CALL condamage(zxx ,zyy ,zzz ,zxy ,zzx ,zyz ,
357 . damage ,nel ,hu ,shape ,whysmax ,
358 . numtabl ,matparam%TABLE ,nvartmp ,vartmp )
359C
360C apply damage
361C
362 DO i=1,nel
363 uvar(i,13) = whysmax(i)
364 DO ik=1,6
365 spk(i,ik)=spknodam(i,ik)*(one-damage(i))
366 ENDDO
367C
368C transform 2PK to Cauchy stresses
369C
370 s1(i) = spk(i,1)
371 s2(i) = spk(i,2)
372 s3(i) = spk(i,3)
373 s4(i) = spk(i,4)
374 s5(i) = spk(i,5)
375 s6(i) = spk(i,6)
376C
377 fs11(i) = f(i,1,1)*s1(i) + f(i,1,2)*s4(i) + f(i,1,3)*s6(i)
378 fs12(i) = f(i,1,1)*s4(i) + f(i,1,2)*s2(i) + f(i,1,3)*s5(i)
379 fs13(i) = f(i,1,1)*s6(i) + f(i,1,2)*s5(i) + f(i,1,3)*s3(i)
380 fs21(i) = f(i,2,1)*s1(i) + f(i,2,2)*s4(i) + f(i,2,3)*s6(i)
381 fs22(i) = f(i,2,1)*s4(i) + f(i,2,2)*s2(i) + f(i,2,3)*s5(i)
382 fs23(i) = f(i,2,1)*s6(i) + f(i,2,2)*s5(i) + f(i,2,3)*s3(i)
383 fs31(i) = f(i,3,1)*s1(i) + f(i,3,2)*s4(i) + f(i,3,3)*s6(i)
384 fs32(i) = f(i,3,1)*s4(i) + f(i,3,2)*s2(i) + f(i,3,3)*s5(i)
385 fs33(i) = f(i,3,1)*s6(i) + f(i,3,2)*s5(i) + f(i,3,3)*s3(i)
386C
387 signxx(i) = fs11(i)*f(i,1,1)+fs12(i)*f(i,1,2)+fs13(i)*f(i,1,3)
388 signyy(i) = fs21(i)*f(i,2,1)+fs22(i)*f(i,2,2)+fs23(i)*f(i,2,3)
389 signzz(i) = fs31(i)*f(i,3,1)+fs32(i)*f(i,3,2)+fs33(i)*f(i,3,3)
390 signxy(i) = fs11(i)*f(i,2,1)+fs12(i)*f(i,2,2)+fs13(i)*f(i,2,3)
391 signyz(i) = fs21(i)*f(i,3,1)+fs22(i)*f(i,3,2)+fs23(i)*f(i,3,3)
392 signzx(i) = fs11(i)*f(i,3,1)+fs12(i)*f(i,3,2)+fs13(i)*f(i,3,3)
393C
394 signxx(i) = signxx(i) / jac(i)
395 signyy(i) = signyy(i) / jac(i)
396 signzz(i) = signzz(i) / jac(i)
397 signxy(i) = signxy(i) / jac(i)
398 signyz(i) = signyz(i) / jac(i)
399 signzx(i) = signzx(i) / jac(i)
400C
401C timestep if the slope of the stress/strain curve
402C exceeds the user specified modulus
403C
404 gs=g
405 ks=k
406 slopemax(i)=max(slope(i,1) , slope(i,2))
407 slopemax(i)=max(slopemax(i), slope(i,3))
408 slopemax(i)=max(slopemax(i), emod)
409
410 gs = slopemax(i)*half
411 ks = slopemax(i)*third
412 !
413 soundsp(i) = sqrt(slopemax(i)/rho(i))
414 viscmax(i) = zero
415 !
416 et(i) = max(slope(i,1),slope(i,2),slope(i,3))/emod
417 ENDDO
418C-------------------------------------------------------------------------
419 RETURN
420 END
max
#define max(a, b)
Definition macros.h:21
condamage_mod
Definition condamage.F:29
condamage_mod::condamage
subroutine condamage(zxx, zyy, zzz, zxy, zzx, zyz, damage, nel, hu, shape, whysmax, numtabl, table, nvartmp, vartmp)
Definition condamage.F:45
conversion_mod
Definition conversion.F:29
conversion_mod::conversion
subroutine conversion(zxx, zyy, zzz, zxy, zyz, zzx, cijkl, dfirst, drate, xscal, scal, nel, jac, slope, numtabl, table, nvartmp, vartmp, yld, ndim_table)
Definition conversion.F:46
prodata
subroutine prodata(a, c, e, nel)
Definition prodATA.F:30
sigeps190
subroutine sigeps190(nel, nuvar, rho, et, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, signxx, signyy, signzz, signxy, signyz, signzx, mfxx, mfxy, mfxz, mfyx, mfyy, mfyz, mfzx, mfzy, mfzz, soundsp, viscmax, uvar, numtabl, matparam, nvartmp, vartmp)
Definition sigeps190.F:45
valpvecdp_v
subroutine valpvecdp_v(sig, val, vec, nel)
Definition sigeps33.F:2231
valpvec_v
subroutine valpvec_v(sig, val, vec, nel)
Definition sigeps33.F:1902