OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps45c.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!|| sigeps45c ../engine/source/materials/mat/mat045/sigeps45c.F
25!||--- called by ------------------------------------------------------
26!|| mulawc ../engine/source/materials/mat_share/mulawc.F90
27!||--- calls -----------------------------------------------------
28!|| finter ../engine/source/tools/curve/finter.f
29!||====================================================================
30 SUBROUTINE sigeps45c(
31 1 NEL0 ,NUPARAM,NUVAR ,NFUNC ,IFUNC ,
32 2 NPF ,NPT0 ,IPT ,IFLAG ,
33 2 TF ,TIME ,TIMESTEP,UPARAM ,RHO0 ,
34 3 AREA ,EINT ,THKLY ,
35 4 EPSPXX ,EPSPYY ,EPSPXY ,EPSPYZ ,EPSPZX ,
36 5 DEPSXX ,DEPSYY ,DEPSXY ,DEPSYZ ,DEPSZX ,
37 6 EPSXX ,EPSYY ,EPSXY ,EPSYZ ,EPSZX ,
38 7 SIGOXX ,SIGOYY ,SIGOXY ,SIGOYZ ,SIGOZX ,
39 8 SIGNXX ,SIGNYY ,SIGNXY ,SIGNYZ ,SIGNZX ,
40 9 SIGVXX ,SIGVYY ,SIGVXY ,SIGVYZ ,SIGVZX ,
41 A SOUNDSP,VISCMAX,THK ,PLA ,UVAR ,
42 B OFF ,NGL ,SHF)
43C-----------------------------------------------
44C I m p l i c i t T y p e s
45C-----------------------------------------------
46#include "implicit_f.inc"
47C-----------------------------------------------
48C G l o b a l P a r a m e t e r s
49C-----------------------------------------------
50#include "mvsiz_p.inc"
51C---------+---------+---+---+--------------------------------------------
52C VAR | SIZE |TYP| RW| DEFINITION
53C---------+---------+---+---+--------------------------------------------
54C NEL0 | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL0
55C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
56C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
57C---------+---------+---+---+--------------------------------------------
58C NFUNC | 1 | I | R | NUMBER FUNCTION USED FOR THIS USER LAW
59C IFUNC | NFUNC | I | R | FUNCTION INDEX
60C NPF | * | I | R | FUNCTION ARRAY
61C NPT0 | 1 | I | R | NUMBER OF LAYERS OR INTEGRATION POINTS
62C IPT | 1 | I | R | LAYER OR INTEGRATION POINT NUMBER
63C IFLAG | * | I | R | GEOMETRICAL FLAGS
64C TF | * | F | R | FUNCTION ARRAY
65C---------+---------+---+---+--------------------------------------------
66C TIME | 1 | F | R | CURRENT TIME
67C TIMESTEP| 1 | F | R | CURRENT TIME STEP
68C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
69C RHO0 | NEL0 | F | R | INITIAL DENSITY
70C AREA | NEL0 | F | R | AREA
71C EINT | 2*NEL0 | F | R | INTERNAL ENERGY(MEMBRANE,BENDING)
72C THKLY | NEL0 | F | R | LAYER THICKNESS
73C EPSPXX | NEL0 | F | R | STRAIN RATE XX
74C EPSPYY | NEL0 | F | R | STRAIN RATE YY
75C ... | | | |
76C DEPSXX | NEL0 | F | R | STRAIN INCREMENT XX
77C DEPSYY | NEL0 | F | R | STRAIN INCREMENT YY
78C ... | | | |
79C EPSXX | NEL0 | F | R | STRAIN XX
80C EPSYY | NEL0 | F | R | STRAIN YY
81C ... | | | |
82C SIGOXX | NEL0 | F | R | OLD ELASTO PLASTIC STRESS XX
83C SIGOYY | NEL0 | F | R | OLD ELASTO PLASTIC STRESS YY
84C ... | | | |
85C---------+---------+---+---+--------------------------------------------
86C SIGNXX | NEL0 | F | W | NEW ELASTO PLASTIC STRESS XX
87C SIGNYY | NEL0 | F | W | NEW ELASTO PLASTIC STRESS YY
88C ... | | | |
89C SIGVXX | NEL0 | F | W | VISCOUS STRESS XX
90C SIGVYY | NEL0 | F | W | VISCOUS STRESS YY
91C ... | | | |
92C SOUNDSP | NEL0 | F | W | SOUND SPEED (NEEDED FOR TIME STEP)
93C VISCMAX | NEL0 | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
94C---------+---------+---+---+--------------------------------------------
95C THK | NEL0 | F |R/W| THICKNESS
96C PLA | NEL0 | F |R/W| PLASTIC STRAIN
97C UVAR |NEL0*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
98C OFF | NEL0 | F |R/W| DELETED ELEMENT FLAG (=1. ON, =0. OFF)
99C---------+---------+---+---+--------------------------------------------
100C-----------------------------------------------
101C I N P U T A r g u m e n t s
102C-----------------------------------------------
103C
104C
105 INTEGER NEL0, NUPARAM, NUVAR, NPT0, IPT,IFLAG(*),
106 . NGL(NEL0)
107 my_real
108 . TIME,TIMESTEP,UPARAM(NUPARAM),
109 . AREA(NEL0),RHO0(NEL0),EINT(NEL0,2),
110 . THKLY(NEL0),PLA(NEL0),
111 . EPSPXX(NEL0),EPSPYY(NEL0),
112 . EPSPXY(NEL0),EPSPYZ(NEL0),EPSPZX(NEL0),
113 . DEPSXX(NEL0),DEPSYY(NEL0),
114 . DEPSXY(NEL0),DEPSYZ(NEL0),DEPSZX(NEL0),
115 . EPSXX(NEL0) ,EPSYY(NEL0) ,
116 . EPSXY(NEL0) ,EPSYZ(NEL0) ,EPSZX(NEL0) ,
117 . sigoxx(nel0),sigoyy(nel0),
118 . sigoxy(nel0),sigoyz(nel0),sigozx(nel0), shf(*)
119C-----------------------------------------------
120C O U T P U T A r g u m e n t s
121C-----------------------------------------------
122 my_real
123 . signxx(nel0),signyy(nel0),
124 . signxy(nel0),signyz(nel0),signzx(nel0),
125 . sigvxx(nel0),sigvyy(nel0),
126 . sigvxy(nel0),sigvyz(nel0),sigvzx(nel0),
127 . soundsp(nel0),viscmax(nel0)
128C-----------------------------------------------
129C I N P U T O U T P U T A r g u m e n t s
130C-----------------------------------------------
131 my_real uvar(nel0,nuvar), off(nel0),thk(nel0)
132C-----------------------------------------------
133C VARIABLES FOR FUNCTION INTERPOLATION
134C-----------------------------------------------
135 INTEGER NPF(*), NFUNC, IFUNC(NFUNC)
136 my_real FINTER ,TF(*)
137 EXTERNAL FINTER
138C Y = FINTER(IFUNC(J),X,NPF,TF,DYDX)
139C Y : y = f(x)
140C X : x
141C DYDX : f'(x) = dy/dx
142C IFUNC(J): FUNCTION INDEX
143C J : FIRST(J=1), SECOND(J=2) .. FUNCTION USED FOR THIS LAW
144C NPF,TF : FUNCTION PARAMETER
145C-----------------------------------------------
146C L o c a l V a r i a b l e s
147C-----------------------------------------------
148 INTEGER I,J,INDEX(MVSIZ),NINDX,NMAX,N
149 my_real
150 . young,anu,g,ca,cb,cn,epsm,sigm,cc,cd,cm,eps0,
151 . ce,ck,c1,c14g3,a1,a2,g3,gs,
152 . ch1,ch2,ch3,qh1,qh2,
153 . nnu1,nu1,nu2,nu3,s1,s2,s3,
154 . r,umr,dezz,
155 . p2,q2,f,df,yld_i,
156 . cutfre,beta
157 my_real
158 . svm(mvsiz),aa(mvsiz),bb(mvsiz),dpla_j(mvsiz),
159 . dpla_i(mvsiz),dr(mvsiz),pp(mvsiz),qq(mvsiz)
160C
161 DATA nmax/3/
162C=======================================================================
163C
164C ZHAO CONSTITUTIVE LAW
165C
166C=======================================================================
167C-----------------------------------------------
168C PARAMETERS READING
169C-----------------------------------------------
170 young = uparam(1)
171 anu = uparam(2)
172 g = uparam(3)
173 ca = uparam(4)
174 cb = uparam(5)
175 cn = uparam(6)
176 epsm = uparam(7)
177 sigm = uparam(8)
178 cc = uparam(9)
179 cd = uparam(10)
180 cm = uparam(11)
181 eps0 = uparam(12)
182 ce = uparam(13)
183 ck = uparam(14)
184 c1 = uparam(15)
185 c14g3 = uparam(16)
186 a1 = uparam(17)
187 a2 = uparam(18)
188 cutfre = uparam(19)
189C-----------------------------------------------
190C USER VARIABLES INITIALIZATION
191C-----------------------------------------------
192 IF(time==zero)THEN
193 DO i=1,nel0
194 uvar(i,1)=zero
195 uvar(i,2)=zero
196 uvar(i,3)=zero
197 uvar(i,4)=zero
198 uvar(i,5)=zero
199 ENDDO
200 ENDIF
201C-----------------------------------------------
202 beta = timestep*2.*pi*cutfre
203 beta = min(one,beta)
204 g3 = three * g
205 gs = five_over_6 * g
206 nnu1 = anu / (one - anu)
207 nu1 = one/(one-anu)
208 nu2 = one/(one+anu)
209 nu3 = one - nnu1
210C NNU2 = NNU1*NNU1
211C NU4 = 1. + NNU2 + NNU1
212C NU5 = 1. + NNU2 - 2.*NNU1
213C NU6 = 0.5- NNU2 + 0.5*NNU1
214C-----------------------------------------------
215C ELASTIC SOLUTION
216C--------------------------------------------
217 DO i=1,nel0
218C
219 signxx(i)=sigoxx(i)+a1*depsxx(i)+a2*depsyy(i)
220 signyy(i)=sigoyy(i)+a2*depsxx(i)+a1*depsyy(i)
221 signxy(i)=sigoxy(i)+g *depsxy(i)
222 signyz(i)=sigoyz(i)+gs *depsyz(i)
223 signzx(i)=sigozx(i)+gs *depszx(i)
224C
225 soundsp(i) = sqrt(a1/rho0(i))
226 viscmax(i) = zero
227C ETSE(I) = 1.
228C-------------------
229C STRAIN RATE
230C-------------------
231C
232 uvar(i,2) = half*( abs(epspxx(i)+epspyy(i))
233 . + sqrt( (epspxx(i)-epspyy(i))*(epspxx(i)-epspyy(i))
234 . + epspxy(i)*epspxy(i) ) )
235 uvar(i,2) = beta*uvar(i,2) + (1.-beta)*uvar(i,5)
236 uvar(i,5) = uvar(i,2)
237C-------------------
238C STRAIN
239C-------------------
240C
241C UVAR(I,1) = 0.5*( EPSXX(I)+EPSYY(I)
242C . + SQRT( (EPSXX(I)-EPSYY(I))*(EPSXX(I)-EPSYY(I))
243C . + EPSXY(I)*EPSXY(I) ) )
244C FAIL(I) = MAX(ZERO,MIN(ONE,(EPSR2(I)-EPST)/(EPSR2(I)-EPSR1(I))))
245C
246 ENDDO
247C-------------------
248C CRITERIA
249C-------------------
250 DO i=1,nel0
251 IF(uvar(i,1)<=zero) THEN
252 ch1=ca
253 ELSEIF(uvar(i,1)>epsm) THEN
254 ch1=ca+cb*epsm**cn
255 ELSE
256 ch1=ca+cb*uvar(i,1)**cn
257 ENDIF
258 IF(uvar(i,2)<=eps0) THEN
259 ch2=0.
260 ELSEIF(uvar(i,1)<=zero) THEN
261 ch2=cc*log(uvar(i,2)/eps0)
262 ELSE
263 ch2=(cc-cd*uvar(i,1)**cm)*log(uvar(i,2)/eps0)
264 ENDIF
265 IF(uvar(i,2)<=zero) THEN
266 ch3=zero
267 ELSE
268 ch3=ce*uvar(i,2)**ck
269 ENDIF
270c jbm033 UVAR(I,3)=MIN(SIGM,CH1+CH2+CH3)
271 uvar(i,3)=min(sigm+ch3,ch1+ch2+ch3)
272 ENDDO
273C------------------------
274C HARDENING MODULUS
275C------------------------
276 DO i=1,nel0
277 IF(uvar(i,1)>zero. and .cn>=one) THEN
278 qh1= cb*cn*uvar(i,1)**(cn-one)
279 ELSEIF(uvar(i,1)>zero. and .cn<one)THEN
280 qh1= cb*cn*uvar(i,1)**(one - cn)
281 ELSE
282 qh1=zero
283 ENDIF
284 IF(uvar(i,1)<=zero. or .uvar(i,2)<=eps0) THEN
285 qh2=zero
286 ELSEIF(cm>=one) THEN
287 qh2=cd*cm*uvar(i,1)**(cm- one)*log(uvar(i,2)/eps0)
288 ELSE
289 qh2=cd*cm*uvar(i,1)**(one -cm)*log(uvar(i,2)/eps0)
290 ENDIF
291 uvar(i,4)=qh1+qh2
292 ENDDO
293C-------------------
294C PROJECTION
295C-------------------
296 IF(iflag(1)==0)THEN
297C radial return
298 DO i=1,nel0
299 svm(i)=sqrt(signxx(i)*signxx(i)
300 . +signyy(i)*signyy(i)
301 . -signxx(i)*signyy(i)
302 . +three*signxy(i)*signxy(i))
303 r = min(one,uvar(i,3)/max(em20,svm(i)))
304 signxx(i)=signxx(i)*r
305 signyy(i)=signyy(i)*r
306 signxy(i)=signxy(i)*r
307 umr = one - r
308 dpla_i(i) = off(i)*svm(i)*umr/(g3+uvar(i,4))
309 uvar(i,1) = uvar(i,1) + dpla_i(i)
310 s1=half*(signxx(i)+signyy(i))
311 dezz = dpla_i(i) * half*(signxx(i)+signyy(i)) /uvar(i,3)
312 dezz=-(depsxx(i)+depsyy(i))*nnu1-nu3*dezz
313 thk(i) = thk(i) + dezz*thkly(i)
314C IF(R<1.) ETSE(I)= H(I)/(H(I)+E(I))
315 ENDDO
316C
317 ELSEIF(iflag(1)==1)THEN
318C-------------------------
319C CRITERE DE VON MISES
320C-------------------------
321 DO i=1,nel0
322 uvar(i,4) = max(zero,uvar(i,4))
323 s1=signxx(i)+signyy(i)
324 s2=signxx(i)-signyy(i)
325 s3=signxy(i)
326 aa(i)=fourth*s1*s1
327 bb(i)=three_over_4*s2*s2+3.*s3*s3
328 svm(i)=sqrt(aa(i)+bb(i))
329 dezz = -(depsxx(i)+depsyy(i))*nnu1
330 thk(i) = thk(i) + dezz*thkly(i)
331 ENDDO
332C-------------------------
333C GATHER PLASTIC FLOW
334C-------------------------
335 nindx=0
336 DO i=1,nel0
337 IF(svm(i)>uvar(i,3).AND.off(i)==one) THEN
338 nindx=nindx+1
339 index(nindx)=i
340 ENDIF
341 ENDDO
342 IF(nindx==0) RETURN
343C---------------------------
344C DEP EN CONTRAINTE PLANE
345C---------------------------
346 DO j=1,nindx
347 i=index(j)
348 dpla_j(i)=(svm(i)-uvar(i,3))/(g3+uvar(i,4))
349C ETSE(I)= H(I)/(H(I)+E(I))
350 ENDDO
351C
352 DO n=1,nmax
353C#include "vectorize.inc"
354 DO j=1,nindx
355 i=index(j)
356 dpla_i(i)=dpla_j(i)
357 yld_i =uvar(i,3)+uvar(i,4)*dpla_i(i)
358 dr(i) =half*young*dpla_i(i)/yld_i
359 pp(i) =one/(one+dr(i)*nu1)
360 qq(i) =one/(one + three*dr(i)*nu2)
361 p2 =pp(i)*pp(i)
362 q2 =qq(i)*qq(i)
363 f =aa(i)*p2+bb(i)*q2-yld_i*yld_i
364 df =-(aa(i)*nu1*p2*pp(i)+three*bb(i)*nu2*q2*qq(i))
365 . *(young-two*dr(i)*uvar(i,4))/yld_i
366 . -two*uvar(i,4)*yld_i
367 IF(dpla_i(i)>zero) THEN
368 dpla_j(i)=max(zero,dpla_i(i)-f/df)
369 ELSE
370 dpla_j(i)=zero
371 ENDIF
372 ENDDO
373 ENDDO
374C------------------------------------------
375C CONTRAINTES PLASTIQUEMENT ADMISSIBLES
376C------------------------------------------
377C#include "vectorize.inc"
378 DO j=1,nindx
379 i=index(j)
380 uvar(i,1) = uvar(i,1) + dpla_i(i)
381 s1=(signxx(i)+signyy(i))*pp(i)
382 s2=(signxx(i)-signyy(i))*qq(i)
383 signxx(i)=half*(s1+s2)
384 signyy(i)=half*(s1-s2)
385 signxy(i)=signxy(i)*qq(i)
386 dezz = - nu3*dr(i)*s1/young
387 thk(i) = thk(i) + dezz*thkly(i)
388 ENDDO
389C
390 ENDIF
391C
392 DO i=1,nel0
393 IF(uvar(i,1)>epsm.AND.off(i)==one)off(i)=four_over_5
394 ENDDO
395C-------------------------
396 RETURN
397 END
398
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
sigeps45c
subroutine sigeps45c(nel0, nuparam, nuvar, nfunc, ifunc, npf, npt0, ipt, iflag, tf, time, timestep, uparam, rho0, area, eint, thkly, epspxx, epspyy, epspxy, epspyz, epspzx, depsxx, depsyy, depsxy, depsyz, depszx, epsxx, epsyy, epsxy, epsyz, epszx, sigoxx, sigoyy, sigoxy, sigoyz, sigozx, signxx, signyy, signxy, signyz, signzx, sigvxx, sigvyy, sigvxy, sigvyz, sigvzx, soundsp, viscmax, thk, pla, uvar, off, ngl, shf)
Definition sigeps45c.F:43