OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps35c.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!|| sigeps35c ../engine/source/materials/mat/mat035/sigeps35c.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 sigeps35c(
31 1 NEL ,NUPARAM,NUVAR ,NFUNC ,IFUNC ,
32 2 NPF ,NPT ,IPT ,IFLA ,
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 ,UVAR ,OFF ,
42 B NGL ,SHF ,ETSE ,ISRATE ,ASRATE ,
43 C EPSD )
44C-----------------------------------------------
45C I m p l i c i t T y p e s
46C-----------------------------------------------
47#include "implicit_f.inc"
48C-----------------------------------------------
49C G l o b a l P a r a m e t e r s
50C-----------------------------------------------
51#include "mvsiz_p.inc"
52C----------------------------------------------------------------
53C I N P U T A R G U M E N T S
54C----------------------------------------------------------------
55 INTEGER NEL,NPT,IPT,IFLA,NUPARAM, NUVAR,NGL(NEL),ISRATE
56
57 my_real
58 . TIME , TIMESTEP , UPARAM(NUPARAM),
59 . RHO0 (NEL), EINT(NEL,2),
60 . EPSPXX(NEL), EPSPYY(NEL),SHF(NEL),
61 . EPSPXY(NEL), EPSPYZ(NEL), EPSPZX(NEL),
62 . DEPSXX(NEL), DEPSYY(NEL),
63 . DEPSXY(NEL), DEPSYZ(NEL), DEPSZX(NEL),
64 . EPSXX (NEL), EPSYY (NEL),
65 . EPSXY (NEL), EPSYZ (NEL), EPSZX (NEL),
66 . SIGOXX(NEL), SIGOYY(NEL),
67 . SIGOXY(NEL), SIGOYZ(NEL), SIGOZX(NEL),
68 . thkly(nel) ,thk(nel),area(nel),
69 . asrate
70C----------------------------------------------------------------
71C O U T P U T A R G U M E N T S
72C----------------------------------------------------------------
73 my_real
74 . signxx(nel), signyy(nel),
75 . signxy(nel), signyz(nel), signzx(nel),
76 . sigvxx(nel), sigvyy(nel),
77 . sigvxy(nel), sigvyz(nel), sigvzx(nel),
78 . soundsp(nel), viscmax(nel),etse(nel)
79C----------------------------------------------------------------
80C I N P U T O U T P U T A R G U M E N T S
81C----------------------------------------------------------------
82 my_real uvar(nel,nuvar), off(nel)
83 my_real ,INTENT(INOUT) :: epsd(nel)
84C----------------------------------------------------------------
85C VARIABLES FOR FUNCTION INTERPOLATION
86C----------------------------------------------------------------
87 INTEGER NPF(*), NFUNC, IFUNC(NFUNC)
88 my_real
89 . FINTER,TF(*)
90 EXTERNAL FINTER
91C----------------------------------------------------------------
92C L O C A L V A R I B L E S
93C----------------------------------------------------------------
94 INTEGER I,KF,IFLAG,ICORRECT
95
96 my_real
97 . poisson,e,e1,e2,g02,
98 . gt2,bulkt,relvexp,
99 . vmu2,vbulk,fac,
100 . c1,c2,c3,pmin,dpdmu,epspc,
101 . p0,phi,gama0, var,epsp,
102 . gama(mvsiz),amu(mvsiz),
103 . sm(mvsiz),em(mvsiz),dedtm(mvsiz),
104 . g2(mvsiz),bulk(mvsiz),
105 . dsxx(mvsiz),dsyy(mvsiz),dsxy(mvsiz),
106 . dexx(mvsiz),deyy(mvsiz),dezz,
107 . dexy(mvsiz),epspzz(mvsiz),epszz(mvsiz),
108 . dedtxx(mvsiz),dedtyy(mvsiz),dedtxy(mvsiz),
109 . dsdtxx(mvsiz),dsdtyy(mvsiz),dsdtxy(mvsiz),
110 . dpdro(mvsiz),p(mvsiz),pdot(mvsiz),
111 . mg2(mvsiz),pg2(mvsiz),mk(mvsiz),pk(mvsiz),
112 . sigair(mvsiz),relvol(mvsiz),enew(mvsiz),
113 . rho(mvsiz)
114
115 my_real
116 . ssm,epsm,epspm,small,cshear,midstep,dt05
117C=======================================================================
118 small = em3
119 cshear= zep426667
120C -- INITIALIZATION---------
121 IF (time == zero) THEN
122 DO i=1,nel
123 uvar(i,3) = area(i)*thk(i) ! VOLUME0
124 ENDDO
125 ENDIF
126C SET INITIAL MATERIAL CONSTANTS
127
128 poisson= uparam(2)
129 e = uparam(1)
130 g02 = half*e/(1+poisson)
131 e1 = uparam(3)
132 e2 = uparam(4)
133 gt2 = uparam(5)/(one+uparam(6))
134 bulkt = uparam(5)/(one-two*uparam(6))
135 vmu2 = two*uparam(7)
136 vbulk = three*uparam(8)+vmu2
137
138 p0 = uparam(9)
139 phi = uparam(10)
140 gama0 = uparam(11)
141
142 c1 = uparam(12)
143 c2 = uparam(13)
144 c3 = uparam(14)
145 iflag = uparam(15)
146 pmin = uparam(16)
147 relvexp = uparam(17)
148 fac = uparam(18)
149
150 kf = ifunc(1)
151C
152 icorrect=0
153 IF(nuparam>=19) icorrect=nint(uparam(19))
154C
155 IF(icorrect == 0)THEN
156 dt05=half*timestep
157 ELSE
158 dt05=zero
159 END IF
160C
161 DO i=1,nel
162 epsxx(i)=epsxx(i)-dt05*epspxx(i)
163 epsyy(i)=epsyy(i)-dt05*epspyy(i)
164 epsxy(i)=epsxy(i)-dt05*epspxy(i)
165 epsyz(i)=epsyz(i)-dt05*epspyz(i)
166 epszx(i)=epszx(i)-dt05*epspzx(i)
167 END DO
168C
169 DO 200 i=1,nel
170c-------calcul E_eq-----
171 epspm = epspxx(i)+epspyy(i)
172 epspc = epspxx(i)-epspyy(i)
173 epsp = sqrt(half*(epspc*epspc+epspxx(i)*epspxx(i)
174 1 + epspyy(i)*epspyy(i)) +three_half* (epspxy(i)*epspxy(i)
175 2 + epspyz(i)*epspyz(i)+epspzx(i)*epspzx(i)))
176 IF (israte > 0) THEN
177 epsd(i) = asrate*epsp + (one - asrate)*uvar(i,4)
178 uvar(i,4) = epsd(i)
179 ELSE
180 epsd(i) = epsp
181 ENDIF
182c
183c------UPDATE ELASTIC MODULUS
184 rho(i)=uvar(i,3)*rho0(i)/area(i)/thk(i)
185 relvol(i)=rho0(i)/rho(i)
186 enew(i)=(max(e,e1*epsp+e2))/(exp(relvexp*log(relvol(i))))
187 g2(i)=enew(i)/(one+poisson)
188 bulk(i)=enew(i)/(one-two*poisson)
189 mk(i) = (bulk(i)+bulkt)/vbulk
190 pk(i) = bulk(i)*bulkt/vbulk
191 mg2(i) = (g2(i)+gt2)/vmu2
192 pg2(i) = g2(i)*gt2/vmu2
193 gama(i) = (relvol(i)-one+gama0)
194 IF(one+gama(i)-phi<=small) gama(i)=-(one-phi-small)
195 ssm=sigoxx(i)+sigoyy(i)
196 sm(i)=third*ssm
197 epsm =epsxx(i)+epsyy(i)
198 var=two*pg2(i)+c3*pk(i)
199 epspzz(i)=((g2(i)-c1*bulk(i))*epspm-(mg2(i)-c2*mk(i))*ssm
200 1 +(pg2(i)-c3*pk(i))*epsm-var*uvar(i,2))
201 2 /(g2(i)+c1*bulk(i)+g2(i)+var*timestep)
202 dezz=epspzz(i)*timestep
203 epszz(i)=uvar(i,2)+dezz
204 uvar(i,2)=epszz(i)
205 thk(i) = thk(i) + dezz*thkly(i)*off(i)
206 em(i)=third*(epsm+epszz(i))
207 dedtm(i)=third*(epspm+epspzz(i))
208 sigair(i)=max(zero,-(p0*gama(i))/(1+gama(i)-phi))
209 200 CONTINUE
210c
211 IF(icorrect/=0)THEN
212 DO 210 i=1,nel
213C COMPUTE DEVIATORIC STRESSES
214 dsxx(i)=sigoxx(i)-sm(i)
215 dsyy(i)=sigoyy(i)-sm(i)
216 dsxy(i)=sigoxy(i)
217
218C COMPUTE DEVIATORIC STRAINS
219 dexx(i)=epsxx(i)-em(i)
220 deyy(i)=epsyy(i)-em(i)
221 dexy(i)=epsxy(i)
222
223C COMPUTE DEVIATORIC STRAIN RATES
224 dedtxx(i)=epspxx(i)-dedtm(i)
225 dedtyy(i)=epspyy(i)-dedtm(i)
226 dedtxy(i)=epspxy(i)
227 210 CONTINUE
228 ELSE
229 DO i=1,nel
230C COMPUTE DEVIATORIC STRESSES
231 dsxx(i)=sigoxx(i)-sm(i)
232 dsyy(i)=sigoyy(i)-sm(i)
233 dsxy(i)=sigoxy(i)
234
235C COMPUTE DEVIATORIC STRAINS
236 dexx(i)=epsxx(i)-em(i)
237 deyy(i)=epsyy(i)-em(i)
238 dexy(i)=epsxy(i)*half
239
240C COMPUTE DEVIATORIC STRAIN RATES
241 dedtxx(i)=epspxx(i)-dedtm(i)
242 dedtyy(i)=epspyy(i)-dedtm(i)
243 dedtxy(i)=epspxy(i)*half
244 END DO
245 END IF
246
247
248C COMPUTE STRESS RATE INCREMENTS
249 DO 250 i=1,nel
250 dsdtxx(i)=g2(i)*dedtxx(i)-mg2(i)*dsxx(i)+pg2(i)*dexx(i)
251 dsdtyy(i)=g2(i)*dedtyy(i)-mg2(i)*dsyy(i)+pg2(i)*deyy(i)
252 dsdtxy(i)=g2(i)*dedtxy(i)-mg2(i)*dsxy(i)+pg2(i)*dexy(i)
253 midstep=one/(one+mg2(i)*dt05)
254 dsdtxx(i)=dsdtxx(i)*midstep
255 dsdtyy(i)=dsdtyy(i)*midstep
256 dsdtxy(i)=dsdtxy(i)*midstep
257 250 CONTINUE
258
259
260 DO 260 i=1,nel
261 sm(i)=sm(i)+uvar(i,1)
262 IF(kf/=0) THEN
263 amu(i)=rho(i)/rho0(i)-one
264 IF(iflag == 0) THEN
265 p(i)=-fac*finter(kf,amu(i),npf,tf,dpdmu)
266 ELSE
267 pdot(i)=( c1*bulk(i)*dedtm(i)
268 & -c2*mk(i)*sm(i)
269 & +c3*pk(i)*em(i))
270 & /(one+c2*dt05*mk(i))
271 p(i)=sm(i)+pdot(i)*timestep
272 & -fac*finter(kf,amu(i),npf,tf,dpdmu)
273 ENDIF
274 ELSE
275 pdot(i)=( c1*bulk(i)*dedtm(i)
276 & -c2*mk(i)*sm(i)
277 & +c3*pk(i)*em(i))
278 & /(one+c2*dt05*mk(i))
279 p(i)=sm(i)+pdot(i)*timestep
280 ENDIF
281 IF(p(i)<=pmin) p(i)=pmin
282 p(i)=p(i)-sigair(i)
283 260 CONTINUE
284
285C COMPUTE SOUND SPEED
286 DO 270 i=1,nel
287 dpdro(i)= g2(i)/(one-poisson)
288 & +p0*(one-phi)/(one+gama(i)-phi)**2
289 270 CONTINUE
290
291 DO 280 i=1,nel
292 soundsp(i)=sqrt(dpdro(i)/rho(i))
293 280 CONTINUE
294C COMPUTE PRESSURE
295
296C COMPUTE UPDATED STRESSES
297 IF(icorrect/=0)THEN
298 DO 290 i=1,nel
299 signxx(i)=dsxx(i)+dsdtxx(i)*timestep+p(i)
300 signyy(i)=dsyy(i)+dsdtyy(i)*timestep+p(i)
301 signxy(i)=dsxy(i)+dsdtxy(i)*timestep*0.5
302 signyz(i)=sigoyz(i)+g02*depsyz(i)*cshear
303 signzx(i)=sigozx(i)+g02*depszx(i)*cshear
304 viscmax(i)= zero
305 uvar(i,1)=sigair(i)
306 290 CONTINUE
307 ELSE
308 DO i=1,nel
309 signxx(i)=dsxx(i)+dsdtxx(i)*timestep+p(i)
310 signyy(i)=dsyy(i)+dsdtyy(i)*timestep+p(i)
311 signxy(i)=dsxy(i)+dsdtxy(i)*timestep
312 signyz(i)=sigoyz(i)+g02*depsyz(i)*cshear
313 signzx(i)=sigozx(i)+g02*depszx(i)*cshear
314 viscmax(i)= zero
315 uvar(i,1)=sigair(i)
316 END DO
317 END IF
318c-----------
319 RETURN
320 END
321
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
max
#define max(a, b)
Definition macros.h:21
sigeps35c
subroutine sigeps35c(nel, nuparam, nuvar, nfunc, ifunc, npf, npt, ipt, ifla, 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, uvar, off, ngl, shf, etse, israte, asrate, epsd)
Definition sigeps35c.F:44