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,J,KF,IFLAG,ICORRECT
95
96 my_real
97 . poisson,e,e1,e2,g02,bulk0,
98 . gt2,bulkt,lamda,lamdat,relvexp,
99 . vmu2,vlamda,vbulk,fac,
100 . c1,c2,c3,pmin,dpdmu,epspc,
101 . p0,phi,gama0, var,epsp,
102 . dpdgama(mvsiz),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,epin,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 epsp = asrate*epsp + (one - asrate)*epsd(i)
178 ENDIF
179 epsd(i) = epsp
180c
181c------UPDATE ELASTIC MODULUS
182 rho(i)=uvar(i,3)*rho0(i)/area(i)/thk(i)
183 relvol(i)=rho0(i)/rho(i)
184 enew(i)=(max(e,e1*epsp+e2))/(exp(relvexp*log(relvol(i))))
185 g2(i)=enew(i)/(one+poisson)
186 bulk(i)=enew(i)/(one-two*poisson)
187 mk(i) = (bulk(i)+bulkt)/vbulk
188 pk(i) = bulk(i)*bulkt/vbulk
189 mg2(i) = (g2(i)+gt2)/vmu2
190 pg2(i) = g2(i)*gt2/vmu2
191 gama(i) = (relvol(i)-one+gama0)
192 IF(one+gama(i)-phi<=small) gama(i)=-(one-phi-small)
193 ssm=sigoxx(i)+sigoyy(i)
194 sm(i)=third*ssm
195 epsm =epsxx(i)+epsyy(i)
196 var=two*pg2(i)+c3*pk(i)
197 epspzz(i)=((g2(i)-c1*bulk(i))*epspm-(mg2(i)-c2*mk(i))*ssm
198 1 +(pg2(i)-c3*pk(i))*epsm-var*uvar(i,2))
199 2 /(g2(i)+c1*bulk(i)+g2(i)+var*timestep)
200 dezz=epspzz(i)*timestep
201 epszz(i)=uvar(i,2)+dezz
202 uvar(i,2)=epszz(i)
203 thk(i) = thk(i) + dezz*thkly(i)*off(i)
204 em(i)=third*(epsm+epszz(i))
205 dedtm(i)=third*(epspm+epspzz(i))
206 sigair(i)=max(zero,-(p0*gama(i))/(1+gama(i)-phi))
207 200 CONTINUE
208c
209 IF(icorrect/=0)THEN
210 DO 210 i=1,nel
211C COMPUTE DEVIATORIC STRESSES
212 dsxx(i)=sigoxx(i)-sm(i)
213 dsyy(i)=sigoyy(i)-sm(i)
214 dsxy(i)=sigoxy(i)
215
216C COMPUTE DEVIATORIC STRAINS
217 dexx(i)=epsxx(i)-em(i)
218 deyy(i)=epsyy(i)-em(i)
219 dexy(i)=epsxy(i)
220
221C COMPUTE DEVIATORIC STRAIN RATES
222 dedtxx(i)=epspxx(i)-dedtm(i)
223 dedtyy(i)=epspyy(i)-dedtm(i)
224 dedtxy(i)=epspxy(i)
225 210 CONTINUE
226 ELSE
227 DO i=1,nel
228C COMPUTE DEVIATORIC STRESSES
229 dsxx(i)=sigoxx(i)-sm(i)
230 dsyy(i)=sigoyy(i)-sm(i)
231 dsxy(i)=sigoxy(i)
232
233C COMPUTE DEVIATORIC STRAINS
234 dexx(i)=epsxx(i)-em(i)
235 deyy(i)=epsyy(i)-em(i)
236 dexy(i)=epsxy(i)*half
237
238C COMPUTE DEVIATORIC STRAIN RATES
239 dedtxx(i)=epspxx(i)-dedtm(i)
240 dedtyy(i)=epspyy(i)-dedtm(i)
241 dedtxy(i)=epspxy(i)*half
242 END DO
243 END IF
244
245
246C COMPUTE STRESS RATE INCREMENTS
247 DO 250 i=1,nel
248 dsdtxx(i)=g2(i)*dedtxx(i)-mg2(i)*dsxx(i)+pg2(i)*dexx(i)
249 dsdtyy(i)=g2(i)*dedtyy(i)-mg2(i)*dsyy(i)+pg2(i)*deyy(i)
250 dsdtxy(i)=g2(i)*dedtxy(i)-mg2(i)*dsxy(i)+pg2(i)*dexy(i)
251 midstep=one/(one+mg2(i)*dt05)
252 dsdtxx(i)=dsdtxx(i)*midstep
253 dsdtyy(i)=dsdtyy(i)*midstep
254 dsdtxy(i)=dsdtxy(i)*midstep
255 250 CONTINUE
256
257
258 DO 260 i=1,nel
259 sm(i)=sm(i)+uvar(i,1)
260 IF(kf/=0) THEN
261 amu(i)=rho(i)/rho0(i)-one
262 IF(iflag == 0) THEN
263 p(i)=-fac*finter(kf,amu(i),npf,tf,dpdmu)
264 ELSE
265 pdot(i)=( c1*bulk(i)*dedtm(i)
266 & -c2*mk(i)*sm(i)
267 & +c3*pk(i)*em(i))
268 & /(one+c2*dt05*mk(i))
269 p(i)=sm(i)+pdot(i)*timestep
270 & -fac*finter(kf,amu(i),npf,tf,dpdmu)
271 ENDIF
272 ELSE
273 pdot(i)=( c1*bulk(i)*dedtm(i)
274 & -c2*mk(i)*sm(i)
275 & +c3*pk(i)*em(i))
276 & /(one+c2*dt05*mk(i))
277 p(i)=sm(i)+pdot(i)*timestep
278 ENDIF
279 IF(p(i)<=pmin) p(i)=pmin
280 p(i)=p(i)-sigair(i)
281 260 CONTINUE
282
283C COMPUTE SOUND SPEED
284 DO 270 i=1,nel
285 dpdro(i)= g2(i)/(one-poisson)
286 & +p0*(one-phi)/(one+gama(i)-phi)**2
287 270 CONTINUE
288
289 DO 280 i=1,nel
290 soundsp(i)=sqrt(dpdro(i)/rho(i))
291 280 CONTINUE
292C COMPUTE PRESSURE
293
294C COMPUTE UPDATED STRESSES
295 IF(icorrect/=0)THEN
296 DO 290 i=1,nel
297 signxx(i)=dsxx(i)+dsdtxx(i)*timestep+p(i)
298 signyy(i)=dsyy(i)+dsdtyy(i)*timestep+p(i)
299 signxy(i)=dsxy(i)+dsdtxy(i)*timestep*0.5
300 signyz(i)=sigoyz(i)+g02*depsyz(i)*cshear
301 signzx(i)=sigozx(i)+g02*depszx(i)*cshear
302 viscmax(i)= zero
303 uvar(i,1)=sigair(i)
304 290 CONTINUE
305 ELSE
306 DO i=1,nel
307 signxx(i)=dsxx(i)+dsdtxx(i)*timestep+p(i)
308 signyy(i)=dsyy(i)+dsdtyy(i)*timestep+p(i)
309 signxy(i)=dsxy(i)+dsdtxy(i)*timestep
310 signyz(i)=sigoyz(i)+g02*depsyz(i)*cshear
311 signzx(i)=sigozx(i)+g02*depszx(i)*cshear
312 viscmax(i)= zero
313 uvar(i,1)=sigair(i)
314 END DO
315 END IF
316c-----------
317 RETURN
318 END
319
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