OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps35.F File Reference
#include "implicit_f.inc"
#include "mvsiz_p.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine sigeps35 (nel, nuparam, nuvar, nfunc, ifunc, npf, tf, time, timestep, uparam, rho0, rho, volume, eint, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, depsxx, depsyy, depszz, depsxy, depsyz, depszx, epsxx, epsyy, epszz, epsxy, epsyz, epszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, sigvxx, sigvyy, sigvzz, sigvxy, sigvyz, sigvzx, soundsp, viscmax, uvar, off, israte, asrate, edot)

Function/Subroutine Documentation

◆ sigeps35()

subroutine sigeps35 ( integer nel,
integer nuparam,
integer nuvar,
integer nfunc,
integer, dimension(nfunc) ifunc,
integer, dimension(*) npf,
tf,
time,
timestep,
uparam,
rho0,
rho,
volume,
eint,
epspxx,
epspyy,
epspzz,
epspxy,
epspyz,
epspzx,
depsxx,
depsyy,
depszz,
depsxy,
depsyz,
depszx,
epsxx,
epsyy,
epszz,
epsxy,
epsyz,
epszx,
sigoxx,
sigoyy,
sigozz,
sigoxy,
sigoyz,
sigozx,
signxx,
signyy,
signzz,
signxy,
signyz,
signzx,
sigvxx,
sigvyy,
sigvzz,
sigvxy,
sigvyz,
sigvzx,
soundsp,
viscmax,
uvar,
off,
integer israte,
asrate,
edot )

Definition at line 31 of file sigeps35.F.

43
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 C o m m o n B l o c k s
54C-----------------------------------------------
55C----------------------------------------------------------------
56C I N P U T A R G U M E N T S
57C----------------------------------------------------------------
58 INTEGER NEL, NUPARAM, NUVAR , ISRATE
59
60 my_real
61 . time , timestep , uparam(nuparam),
62 . rho(nel), rho0(nel), volume(nel), eint(nel),
63 . epspxx(nel), epspyy(nel), epspzz(nel),
64 . epspxy(nel), epspyz(nel), epspzx(nel),
65 . depsxx(nel), depsyy(nel), depszz(nel),
66 . depsxy(nel), depsyz(nel), depszx(nel),
67 . epsxx(nel), epsyy(nel), epszz(nel),
68 . epsxy(nel), epsyz(nel), epszx(nel),
69 . sigoxx(nel), sigoyy(nel), sigozz(nel),
70 . sigoxy(nel), sigoyz(nel), sigozx(nel), asrate
71C----------------------------------------------------------------
72C O U T P U T A R G U M E N T S
73C----------------------------------------------------------------
74 my_real
75 . signxx(nel), signyy(nel), signzz(nel),
76 . signxy(nel), signyz(nel), signzx(nel),
77 . sigvxx(nel), sigvyy(nel), sigvzz(nel),
78 . sigvxy(nel), sigvyz(nel), sigvzx(nel),
79 . soundsp(nel), viscmax(nel)
80C----------------------------------------------------------------
81C I N P U T O U T P U T A R G U M E N T S
82C----------------------------------------------------------------
83 my_real uvar(nel,nuvar), off(nel) , edot(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 . e,poisson,a,b,et,poissont,fac,epsp
98 my_real
99 . gt2,bulkt3,relvexp
100 my_real
101 . vmu,vlamda,vmu2,vlamda3
102 my_real
103 . c1,c2,c3,pmin,dpdmu(mvsiz)
104 my_real
105 . p0,phi,gama0
106 my_real
107 . enew(mvsiz),dedot(mvsiz)
108 my_real
109 . dpdgama(mvsiz),gama(mvsiz),amu(mvsiz)
110 my_real
111 . sm(mvsiz),em(mvsiz),dedtm(mvsiz)
112 my_real
113 . g2(mvsiz),bulk(mvsiz),bulk3(mvsiz)
114 my_real
115 . dsxx(mvsiz),dsyy(mvsiz),dszz(mvsiz)
116 my_real
117 . dsxy(mvsiz),dsyz(mvsiz),dszx(mvsiz)
118 my_real
119 . dexx(mvsiz),deyy(mvsiz),dezz(mvsiz)
120 my_real
121 . dexy(mvsiz),deyz(mvsiz),dezx(mvsiz)
122 my_real
123 . dedtxx(mvsiz),dedtyy(mvsiz),dedtzz(mvsiz)
124 my_real
125 . dedtxy(mvsiz),dedtyz(mvsiz),dedtzx(mvsiz)
126 my_real
127 . dsdtxx(mvsiz),dsdtyy(mvsiz),dsdtzz(mvsiz)
128 my_real
129 . dsdtxy(mvsiz),dsdtyz(mvsiz),dsdtzx(mvsiz)
130 my_real
131 . dpdro(mvsiz),p(mvsiz),pdot(mvsiz),relvol(mvsiz)
132 my_real
133 . sigair(mvsiz)
134
135 my_real
136 . small,tiny , aux1, aux2,
137 . midstep, dt05
138 LOGICAL TEST
139C----------------------------------------------------------------
140 small = em3
141 tiny = em30
142 test=.false.
143
144C SET INITIAL MATERIAL CONSTANTS
145
146 e = uparam(1)
147 poisson = uparam(2)
148 a = uparam(3)
149 b = uparam(4)
150
151 et = uparam(5)
152 poissont= uparam(6)
153 vmu2 = uparam(7)*two
154 vlamda3 = uparam(8)*three
155
156 p0 = uparam(9)
157 phi = uparam(10)
158 gama0 = uparam(11)
159
160 c1 = uparam(12)
161 c2 = uparam(13)
162 c3 = uparam(14)
163
164 iflag = uparam(15)
165 pmin = uparam(16)
166 relvexp = uparam(17)
167
168 fac = uparam(18)
169
170 kf = ifunc(1)
171
172 gt2 = et/(one+poissont)
173 bulkt3 = et/(one-two*poissont)
174
175 icorrect=0
176 IF(nuparam>=19) icorrect=nint(uparam(19))
177C
178 IF(icorrect==0)THEN
179 dt05=half*timestep
180 ELSE
181 dt05=zero
182 END IF
183C
184 DO i=1,nel
185 epsxx(i)=epsxx(i)-dt05*epspxx(i)
186 epsyy(i)=epsyy(i)-dt05*epspyy(i)
187 epszz(i)=epszz(i)-dt05*epspzz(i)
188 epsxy(i)=epsxy(i)-dt05*epspxy(i)
189 epsyz(i)=epsyz(i)-dt05*epspyz(i)
190 epszx(i)=epszx(i)-dt05*epspzx(i)
191 END DO
192
193 DO 200 i=1,nel
194 gama(i) = (rho0(i)/rho(i)-one+gama0)
195 IF(one+gama(i)-phi<=small) gama(i)=-(one-phi-small)
196 sm(i)=third*(sigoxx(i)+sigoyy(i)+sigozz(i))+uvar(i,1)
197 em(i)=third*(epsxx(i)+epsyy(i)+epszz(i))
198 dedtm(i)=third*(epspxx(i)+epspyy(i)+epspzz(i))
199 sigair(i)=max(zero,-(p0*gama(i))/(one+gama(i)-phi))
200 200 CONTINUE
201
202 DO 210 i=1,nel
203C COMPUTE DEVIATORIC STRESSES
204 dsxx(i)=sigoxx(i)-sm(i)+uvar(i,1)
205 dsyy(i)=sigoyy(i)-sm(i)+uvar(i,1)
206 dszz(i)=sigozz(i)-sm(i)+uvar(i,1)
207 dsxy(i)=sigoxy(i)
208 dsyz(i)=sigoyz(i)
209 dszx(i)=sigozx(i)
210
211C COMPUTE DEVIATORIC STRAINS
212 dexx(i)=epsxx(i)-em(i)
213 deyy(i)=epsyy(i)-em(i)
214 dezz(i)=epszz(i)-em(i)
215 dexy(i)=epsxy(i)* half
216 deyz(i)=epsyz(i)* half
217 dezx(i)=epszx(i)* half
218
219C COMPUTE DEVIATORIC STRAIN RATES
220 dedtxx(i)=epspxx(i)-dedtm(i)
221 dedtyy(i)=epspyy(i)-dedtm(i)
222 dedtzz(i)=epspzz(i)-dedtm(i)
223 dedtxy(i)=epspxy(i)*half
224 dedtyz(i)=epspyz(i)*half
225 dedtzx(i)=epspzx(i)*half
226 210 CONTINUE
227
228 DO i=1,nel
229 relvol(i)=rho0(i)/rho(i)
230 epsp = max(
231 & abs(epspxx(i)),abs(epspyy(i)),abs(epspzz(i)),
232 & abs(epspxy(i)),abs(epspyz(i)),abs(epspzx(i)))
233 IF (israte == 0) THEN
234 edot(i) = epsp
235 ELSE
236 edot(i) = asrate*epsp + (one - asrate)*uvar(i,4)
237 uvar(i,4) = edot(i)
238 ENDIF
239 ENDDO
240
241C UPDATE ELASTIC MODULUS
242 IF(relvexp/=zero) THEN
243 DO i=1,nel
244 enew(i)=(max(e,a*edot(i)+b))/(exp(relvexp*log(relvol(i))))
245 ENDDO
246 ELSE
247 DO i=1,nel
248 enew(i)=(max(e,a*edot(i)+b))
249 ENDDO
250 ENDIF
251
252
253 aux1 = one / (one+poisson)
254 aux2 = one / (one-two*poisson)
255 DO i=1,nel
256 g2(i)=enew(i)*aux1
257 bulk3(i)=enew(i)*aux2
258 ENDDO
259C COMPUTE DEVIATORIC STRESS RATE INCREMENTS
260 DO i=1,nel
261 dsdtxx(i)=g2(i)*dedtxx(i)-(g2(i)+gt2)*dsxx(i)/vmu2+
262 & g2(i)*gt2*dexx(i)/vmu2
263 dsdtyy(i)=g2(i)*dedtyy(i)-(g2(i)+gt2)*dsyy(i)/vmu2+
264 & g2(i)*gt2*deyy(i)/vmu2
265 dsdtzz(i)=g2(i)*dedtzz(i)-(g2(i)+gt2)*dszz(i)/vmu2+
266 & g2(i)*gt2*dezz(i)/vmu2
267 dsdtxy(i)=g2(i)*dedtxy(i)-(g2(i)+gt2)*dsxy(i)/vmu2+
268 & g2(i)*gt2*dexy(i)/vmu2
269 dsdtyz(i)=g2(i)*dedtyz(i)-(g2(i)+gt2)*dsyz(i)/vmu2+
270 & g2(i)*gt2*deyz(i)/vmu2
271 dsdtzx(i)=g2(i)*dedtzx(i)-(g2(i)+gt2)*dszx(i)/vmu2+
272 & g2(i)*gt2*dezx(i)/vmu2
273 midstep=one/(one+(g2(i)+gt2)/vmu2*dt05)
274 dsdtxx(i)=dsdtxx(i)*midstep
275 dsdtyy(i)=dsdtyy(i)*midstep
276 dsdtzz(i)=dsdtzz(i)*midstep
277 dsdtxy(i)=dsdtxy(i)*midstep
278 dsdtyz(i)=dsdtyz(i)*midstep
279 dsdtzx(i)=dsdtzx(i)*midstep
280 ENDDO
281
282
283C COMPUTE PRESSURE
284 DO i=1,nel
285 IF(kf/=0) THEN
286 amu(i)=rho(i)/rho0(i)-one
287 IF(iflag==0) THEN
288 p(i)=-fac*finter(kf,amu(i),npf,tf,dpdmu(i))
289 ELSE
290 pdot(i)=( c1*(bulk3(i)*dedtm(i))
291 & -c2*((bulk3(i)+bulkt3)*sm(i)/(vlamda3+vmu2))
292 & +c3*((bulk3(i)*bulkt3)*em(i)/(vlamda3+vmu2)))
293 & / (one+c2*(bulk3(i)+bulkt3)/(vlamda3+vmu2)*dt05)
294 p(i)=sm(i)+pdot(i)*timestep
295 & -fac*finter(kf,amu(i),npf,tf,dpdmu(i))
296 ENDIF
297 ELSE
298 pdot(i)=( c1*(bulk3(i)*dedtm(i))
299 & -c2*((bulk3(i)+bulkt3)*sm(i)/(vlamda3+vmu2))
300 & +c3*((bulk3(i)*bulkt3)*em(i)/(vlamda3+vmu2)))
301 & / (one+c2*(bulk3(i)+bulkt3)/(vlamda3+vmu2)*dt05)
302 p(i)=sm(i)+pdot(i)*timestep
303 ENDIF
304 IF(p(i)<=pmin) p(i)=pmin
305 ENDDO
306
307C COMPUTE SOUND SPEED
308 DO i=1,nel
309 IF(kf==0)THEN
310 dpdro(i)=two_third*g2(i)+third*bulk3(i)
311 & +p0*(one-phi)/(one+gama(i)-phi)**2
312 ELSEIF(iflag==0)THEN
313C
314C BULK may be unnecessary ; stability over 2 cycles is not fully insured.
315C---normally BULK3 is not used in this case, just keep it as assurance
316 dpdro(i)=two_third*g2(i)
317 & +max(third*bulk3(i),abs(fac*dpdmu(i)))
318 & +p0*(one-phi)/(one+gama(i)-phi)**2
319 ELSE
320C
321C stability over 2 cycles is not fully insured.
322 dpdro(i)=two_third*g2(i)+third*bulk3(i)+abs(fac*dpdmu(i))
323 & +p0*(one-phi)/(one+gama(i)-phi)**2
324 END IF
325 ENDDO
326
327 DO i=1,nel
328 soundsp(i)=sqrt(dpdro(i)/rho(i))
329 ENDDO
330
331C COMPUTE UPDATED STRESSES
332 DO i=1,nel
333 signxx(i)=dsxx(i)+dsdtxx(i)*timestep+p(i)-sigair(i)
334 signyy(i)=dsyy(i)+dsdtyy(i)*timestep+p(i)-sigair(i)
335 signzz(i)=dszz(i)+dsdtzz(i)*timestep+p(i)-sigair(i)
336 signxy(i)=dsxy(i)+dsdtxy(i)*timestep
337 signyz(i)=dsyz(i)+dsdtyz(i)*timestep
338 signzx(i)=dszx(i)+dsdtzx(i)*timestep
339 viscmax(i)= zero
340 uvar(i,1)=sigair(i)
341 ENDDO
342
343 RETURN
344
my_real
#define my_real
Definition cppsort.cpp:32
max
#define max(a, b)
Definition macros.h:21