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

Go to the source code of this file.

Functions/Subroutines

subroutine sigeps78 (nel, nuparam, nuvar, nfunc, ifunc, npf, tf, time, timestep, uparam, rho0, rho, siga, sigb, sigc, uvar, pla, dep, depsxx, depsyy, depszz, depsxy, depsyz, depszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, soundsp, off, yld, etse)

Function/Subroutine Documentation

◆ sigeps78()

subroutine sigeps78 ( integer nel,
integer nuparam,
integer nuvar,
integer nfunc,
integer, dimension(nfunc) ifunc,
integer, dimension(*) npf,
tf,
time,
timestep,
uparam,
rho0,
rho,
intent(inout) siga,
intent(inout) sigb,
intent(inout) sigc,
uvar,
pla,
dep,
depsxx,
depsyy,
depszz,
depsxy,
depsyz,
depszx,
sigoxx,
sigoyy,
sigozz,
sigoxy,
sigoyz,
sigozx,
signxx,
signyy,
signzz,
signxy,
signyz,
signzx,
soundsp,
off,
yld,
etse )

Definition at line 33 of file sigeps78.F.

41C-----------------------------------------------
42C I m p l i c i t T y p e s
43C-----------------------------------------------
44#include "implicit_f.inc"
45C-----------------------------------------------
46C C O M M O N
47C-----------------------------------------------
48#include "com01_c.inc"
49C-----------------------------------------------
50C I N P U T A r g u m e n t s
51C-----------------------------------------------
52 my_real ,DIMENSION(6,NEL), INTENT(INOUT) :: siga,sigb,sigc
53 INTEGER NEL, NUPARAM, NUVAR
54 my_real
55 . time,timestep,uparam(nuparam),
56 . rho(nel),rho0(nel),
57 . depsxx(nel),depsyy(nel),depszz(nel),
58 . depsxy(nel),depsyz(nel),depszx(nel),
59 . sigoxx(nel),sigoyy(nel),sigozz(nel),
60 . sigoxy(nel),sigoyz(nel),sigozx(nel)
61C-----------------------------------------------
62C O U T P U T A r g u m e n t s
63C-----------------------------------------------
64 my_real
65 . signxx(nel),signyy(nel),signzz(nel),
66 . signxy(nel),signyz(nel),signzx(nel),
67 . soundsp(nel),yld(nel) ,
68 . pla(nel),dep(nel),etse(nel)
69C-----------------------------------------------
70C I N P U T O U T P U T A r g u m e n t s
71C-----------------------------------------------
72 my_real
73 . uvar(nel,nuvar), off(nel)
74C-----------------------------------------------
75C VARIABLES FOR FUNCTION INTERPOLATION
76C-----------------------------------------------
77 INTEGER NPF(*), NFUNC, IFUNC(NFUNC)
78 my_real
79 . finter ,tf(*)
80 EXTERNAL finter
81C Y = FINTER(IFUNC(J),X,NPF,TF,DYDX)
82C Y : y = f(x)
83C X : x
84C DYDX : f'(x) = dy/dx
85C IFUNC(J): FUNCTION INDEX
86C J : FIRST(J=1), SECOND(J=2) .. FUNCTION USED FOR THIS LAW
87C NPF,TF : FUNCTION PARAMETER
88C-----------------------------------------------
89C L o c a l V a r i a b l e s
90C-----------------------------------------------
91 INTEGER I,J,K,NITER,OPTE
92 my_real eini,nu,bsat,hyu,cyu,rsat,
93 . yield,byu,myu,c1_kh,c2_kh,
94 . einf,coe,fctp,fct,cun,cdeux,ctrois,sigy,dmu,prden,
95 . tha,thb,tolr,seff,dg,seffp,casta,thetaq,dtca,dtbm,
96 . sthabxx,sthabyy,sthabzz,sthabxy,sthabyz,sthabzx,
97 . saxx,sayy,sazz,saxy,sayz,sazx,gama,hdgsig,h,bqdb,eqbq
98 my_real :: young(nel),
99 . alxx(nel),alyy(nel),alzz(nel),rbulk(nel),shear(nel),lamda(nel),
100 . alxy(nel),alyz(nel),alzx(nel),
101 . ayu(nel),sigm(nel),sigeq(nel),dydxe(nel),wave(nel),r(nel),
102 . strialxx(nel),strialyy(nel),strialzz(nel),
103 . strialxy(nel),strialyz(nel),strialzx(nel),
104 . g2(nel),asta(nel),scalee(nel),rnih(nel),drnih(nel),epbar(nel),
105 . max_asta(nel)
106 my_real :: db(6,nel),ep(6,nel),ast(6,nel),beta(6,nel),
107 . dr(nel),bq(6,nel)
108C=======================================================================
109C SIGA = CENTER OF YIELD SURFACE ALPHA
110C SIGB = CENTER OF BOUNDING SURFACE
111C SIGC = q : CENTER OF NON-IH SURFACE G_SIGMA
112c
113C UVAR(I,1) = R : ISOTROPIC HARDENING OF BOUNDING SURFACE
114C UVAR(I,2) = r : RADIUS OF G_SIGMA (RNIH)
115C UVAR(I,3) = a = B + R - YIELD (AYU)
116C UVAR(I,4) = YLD STRESS
117C UVAR(I,5) = EFFECTIVE PLASTIC STRAIN INCREMENT
118C=======================================================================
119 eini = uparam(1)
120 nu = uparam(2)
121 yield = uparam(3)
122 byu = uparam(4)
123 c2_kh = uparam(5)
124 hyu = uparam(6)
125 bsat = uparam(7)
126 myu = uparam(8)
127 rsat = uparam(9)
128 einf = uparam(10)
129 coe = uparam(11)
130 opte = uparam(12)
131 c1_kh = uparam(22)
132c
133 niter = 3
134 cun = third/(one-two*nu)
135 cdeux = half/(one+nu)
136 ctrois = nu/(one+nu)/(one-two*nu)
137 IF (isigi == 0) THEN
138 IF(time == zero) THEN
139 DO i=1,nel
140 uvar(i,3) = bsat - yield
141 uvar(i,4) = yield
142 ENDDO
143 ENDIF
144 ENDIF
145C EINI = initial young modulus
146C EINF = saturated young modulus
147C PLA=equivalent plastic strain (eps bar)
148C READING STORED STATE_VARIABLES
149 DO j=1,6
150 DO i=1,nel
151 ast(j,i) =siga(j,i)
152 beta(j,i)=sigb(j,i)
153 ENDDO
154 ENDDO
155 IF (opte == 1)THEN
156 DO i=1,nel
157 young(i)=eini
158 IF(pla(i) > zero)THEN
159 scalee(i) = finter(ifunc(1),pla(i),npf,tf,dydxe(i))
160 young(i) = scalee(i)*eini
161 ENDIF
162 ENDDO
163 ELSE
164 DO i=1,nel
165 young(i)=eini
166 IF(pla(i) > zero)THEN
167 young(i)=eini-(eini-einf)*(one-exp(-coe*pla(i)))
168 ENDIF
169 ENDDO
170 ENDIF
171 DO i=1,nel
172 rbulk(i) = young(i)*cun
173 shear(i) = young(i)*cdeux
174 lamda(i) = young(i)*ctrois
175 g2(i) = two*shear(i)
176C Module d'onde de compression
177 wave(i)=g2(i)+lamda(i)
178 asta(i)=zero
179 max_asta(i) = uvar(i,6)
180 ENDDO
181
182 DO i=1,nel
183 alxx(i)=ast(1,i)+beta(1,i)
184 alyy(i)=ast(2,i)+beta(2,i)
185 alzz(i)=ast(3,i)+beta(3,i)
186 alxy(i)=ast(4,i)+beta(4,i)
187 alyz(i)=ast(5,i)+beta(5,i)
188 alzx(i)=ast(6,i)+beta(6,i)
189 asta(i)=asta(i)+three_half*(
190 . ast(1,i)*ast(1,i)+ast(2,i)*ast(2,i)+ast(3,i)*ast(3,i)
191 . +two*(ast(4,i)*ast(4,i)+ast(5,i)*ast(5,i)
192 . +ast(6,i)*ast(6,i)))
193 asta(i)=sqrt(max(em20,asta(i)))
194 max_asta(i) = max(max_asta(i),asta(i))
195 ENDDO
196C =====================================================
197C DEVIATORIC STRESS AT CURRENT STEP
198C =====================================================
199 DO i=1,nel
200C Estimation sigma_n+1
201 signxx(i)=sigoxx(i)+wave(i)*depsxx(i)
202 . +lamda(i)*depsyy(i)
203 . +lamda(i)*depszz(i)
204 signyy(i)=sigoyy(i)+wave(i)*depsyy(i)
205 . +lamda(i)*depsxx(i)
206 . +lamda(i)*depszz(i)
207 signzz(i)=sigozz(i)+wave(i)*depszz(i)
208 . +lamda(i)*depsxx(i)
209 . +lamda(i)*depsyy(i)
210 sigm(i)=-(signxx(i)+signyy(i)+signzz(i)) * third
211C STRIAL_ij deviateur de sigma_n+1=sigm_n+ELa DEPS
212 strialxx(i)=signxx(i)+sigm(i)
213 strialyy(i)=signyy(i)+sigm(i)
214 strialzz(i)=signzz(i)+sigm(i)
215 strialxy(i)=sigoxy(i)+g2(i)*depsxy(i)
216 strialyz(i)=sigoyz(i)+g2(i)*depsyz(i)
217 strialzx(i)=sigozx(i)+g2(i)*depszx(i)
218 ENDDO
219C =====================================================
220C CURRENT TRIAL EQUIVALENT STRESS CALCULATION
221C =====================================================
222 DO i=1,nel
223C f=3/2(s-alpha):(s-alpha)-Y.Y
224 seff=sqrt(three_half*(
225 . (strialxx(i)-alxx(i))*(strialxx(i)-alxx(i))
226 . +(strialyy(i)-alyy(i))*(strialyy(i)-alyy(i))
227 . +(strialzz(i)-alzz(i))*(strialzz(i)-alzz(i))
228 . +two*((strialxy(i)-alxy(i))*(strialxy(i)-alxy(i))
229 . +(strialyz(i)-alyz(i))*(strialyz(i)-alyz(i))
230 . +(strialzx(i)-alzx(i))*(strialzx(i)-alzx(i)))))
231 dep(i)=zero
232 ayu(i)=uvar(i,3)
233 IF (seff <= yield) THEN
234 DO j=1,6
235 ep(j,i)=zero
236 ENDDO
237 ELSE
238C=======================================================================
239C------------------------------------------------------------
240C PROJECTION
241C------------------------------------------------------------
242C=======================================================================
243C Compute the effective strain increment with Newton algorithm
244 dep(i)=uvar(i,5)
245C=======================================================================
246C Newton
247C=======================================================================
248 IF (max_asta(i) < bsat - yield) THEN
249 cyu = c1_kh
250 ELSE
251 cyu = c2_kh
252 ENDIF
253 prden=three*shear(i)+cyu*ayu(i)+myu*byu
254 casta=cyu*sqrt(ayu(i)/asta(i))
255 DO k=1,niter
256c THA et THB pour semi implicite---------------------------------
257 tha=one-casta*dep(i)
258 thb=one-myu*dep(i)
259C sum of Tetaa alpha* + Tetab Beta=STHAB
260 sthabxx=tha*ast(1,i)+thb*beta(1,i)
261 sthabyy=tha*ast(2,i)+thb*beta(2,i)
262 sthabzz=tha*ast(3,i)+thb*beta(3,i)
263 sthabxy=tha*ast(4,i)+thb*beta(4,i)
264 sthabyz=tha*ast(5,i)+thb*beta(5,i)
265 sthabzx=tha*ast(6,i)+thb*beta(6,i)
266 seff=sqrt(three_half*(
267 . (strialxx(i)-sthabxx)*(strialxx(i)-sthabxx)
268 . +(strialyy(i)-sthabyy)*(strialyy(i)-sthabyy)
269 . +(strialzz(i)-sthabzz)*(strialzz(i)-sthabzz)
270 . +two*((strialxy(i)-sthabxy)*(strialxy(i)-sthabxy)
271 . +(strialyz(i)-sthabyz)*(strialyz(i)-sthabyz)
272 . +(strialzx(i)-sthabzx)*(strialzx(i)-sthabzx))))
273 seff=max(em20,seff)
274 fct=seff-yield-prden*dep(i)
275 seffp=three*((strialxx(i)-sthabxx)
276 . *(casta*ast(1,i)+myu*beta(1,i))
277 . +(strialyy(i)-sthabyy)
278 . *(casta*ast(2,i)+myu*beta(2,i))
279 . +(strialzz(i)-sthabzz)
280 . *(casta*ast(3,i)+myu*beta(3,i))
281 . +two*((strialxy(i)-sthabxy)
282 . *(casta*ast(4,i)+myu*beta(4,i))
283 . +(strialyz(i)-sthabyz)
284 . *(casta*ast(5,i)+myu*beta(5,i))
285 . +(strialzx(i)-sthabzx)
286 . *(casta*ast(6,i)+myu*beta(6,i))))
287c
288 fctp=half*seffp/seff-prden
289 dep(i)=max(zero,dep(i)-fct/fctp)
290 ENDDO
291C=======================================================================
292C END of Newton iteration
293C=======================================================================
294 tha=one-casta*dep(i)
295 thb=one-myu*dep(i)
296C somme de Tetaa alpha* + Tetab Beta=STHAB
297 sthabxx=tha*ast(1,i)+thb*beta(1,i)
298 sthabyy=tha*ast(2,i)+thb*beta(2,i)
299 sthabzz=tha*ast(3,i)+thb*beta(3,i)
300 sthabxy=tha*ast(4,i)+thb*beta(4,i)
301 sthabyz=tha*ast(5,i)+thb*beta(5,i)
302 sthabzx=tha*ast(6,i)+thb*beta(6,i)
303c
304c SEFF=SQRT(THREE*HALF*(
305c . (STRIALXX(I)-STHABXX)*(STRIALXX(I)-STHABXX)
306c . +(STRIALYY(I)-STHABYY)*(STRIALYY(I)-STHABYY)
307c . +(STRIALZZ(I)-STHABZZ)*(STRIALZZ(I)-STHABZZ)
308c . +TWO*((STRIALXY(I)-STHABXY)*(STRIALXY(I)-STHABXY)
309c . +(STRIALYZ(I)-STHABYZ)*(STRIALYZ(I)-STHABYZ)
310c . +(STRIALZX(I)-STHABZX)*(STRIALZX(I)-STHABZX))))
311c SEFF=MAX(EM20,SEFF)
312c ========================================================================
313C Calcul de s_n+1 - alpha_n+1
314 saxx=yield*(strialxx(i)-sthabxx)/(yield+dep(i)*prden)
315 sayy=yield*(strialyy(i)-sthabyy)/(yield+dep(i)*prden)
316 sazz=yield*(strialzz(i)-sthabzz)/(yield+dep(i)*prden)
317 saxy=yield*(strialxy(i)-sthabxy)/(yield+dep(i)*prden)
318 sayz=yield*(strialyz(i)-sthabyz)/(yield+dep(i)*prden)
319 sazx=yield*(strialzx(i)-sthabzx)/(yield+dep(i)*prden)
320C============= EP= delta epsilon_n+1^p=======================
321 ep(1,i)=three_half*dep(i)*saxx/yield
322 ep(2,i)=three_half*dep(i)*sayy/yield
323 ep(3,i)=three_half*dep(i)*sazz/yield
324 ep(4,i)=three_half*dep(i)*saxy/yield
325 ep(5,i)=three_half*dep(i)*sayz/yield
326 ep(6,i)=three_half*dep(i)*sazx/yield
327c EPBAR(I)=SQRT(2./3.*(EP(1,I)*EP(1,I)+EP(2,I)*EP(2,I)
328c . +EP(3,I)*EP(3,I)+TWO*(EP(4,I)*EP(4,I)+EP(5,I)*EP(5,I)
329c . +EP(6,I)*EP(6,I))))
330C=======================================================================
331 dtca=two*third*cyu*ayu(i)
332 dtbm=two*third*byu*myu
333 DO j=1,6
334C Actualisation de alpha* et beta avec Epsilon_n+1^p
335 ast(j,i)=tha*ast(j,i)+dtca*ep(j,i)
336 beta(j,i)=thb*beta(j,i)+dtbm*ep(j,i)
337 ENDDO
338 DO j=1,6
339C Compute betadot
340 db(j,i)=beta(j,i)-sigb(j,i)
341 siga(j,i)=ast(j,i)
342 sigb(j,i)=beta(j,i)
343 ENDDO
344C=======================================================================
345C WORKHARDENING STAGNATION
346C=======================================================================
347 r(i)=uvar(i,1)
348 bqdb=zero
349 eqbq=zero
350 rnih(i)=uvar(i,2)
351 DO j=1,6
352C BQ= beta_n+1-q_n
353 bq(j,i)=beta(j,i)-sigc(j,i)
354C Compute g_sigma the non-IH surface
355 eqbq=eqbq+three_half*bq(j,i)*bq(j,i)
356 bqdb=bqdb+bq(j,i)*db(j,i)
357 ENDDO
358 tolr=eqbq-rnih(i)*rnih(i)
359 IF (tolr >= zero .AND. bqdb > zero) THEN
360 r(i)=thb*uvar(i,1)+myu*rsat*dep(i) ! Rdot=m(R_sat-R).pdot
361 dr(i)=r(i)-uvar(i,1)
362 ELSE
363 dr(i)=zero
364 ENDIF
365 gama=zero
366 dmu=zero
367 IF ( hyu> zero )THEN
368 hdgsig=three*hyu*bqdb
369 IF (rnih(i) == zero)THEN
370 dmu = eqbq/hdgsig-one
371 ELSE
372 dmu =((three*hyu*bqdb+sqrt(hdgsig*hdgsig
373 . +four*rnih(i)*rnih(i)*eqbq))
374 . /(two*rnih(i)*rnih(i)))-one
375 ENDIF
376 bqdb=zero
377 DO j=1,6
378 sigc(j,i)=beta(j,i)-bq(j,i)/(one+dmu) !update of q_n+1
379 ENDDO
380 DO j=1,6
381 bq(j,i)=beta(j,i)-sigc(j,i)
382 ENDDO
383 DO j=1,6
384 bqdb=bqdb+bq(j,i)*db(j,i)
385 ENDDO
386 IF (dr(i) > zero) THEN
387 IF ( rnih(i) == zero) THEN
388 rnih(i)=three*bqdb*hyu
389 ELSE
390 gama=bqdb*three_half/rnih(i)
391 rnih(i)= uvar(i,2)+ hyu*gama
392 ENDIF
393 ENDIF
394 ENDIF
395 uvar(i,1)=r(i)
396 uvar(i,2)=rnih(i)
397 ayu(i)=bsat+r(i)-yield
398 uvar(i,3)=ayu(i)
399C END WORKHARDENING STAGNATION
400 ENDIF
401 ENDDO
402C=======================================================================
403C UPDATE CURRENT STRESSES
404C============================================================
405 DO i=1,nel
406 signxx(i)=sigoxx(i)+wave(i)*(depsxx(i)-ep(1,i))
407 . +lamda(i)*(depsyy(i)-ep(2,i))
408 . +lamda(i)*(depszz(i)-ep(3,i))
409 signyy(i)=sigoyy(i)+wave(i)*(depsyy(i)-ep(2,i))
410 . +lamda(i)*(depsxx(i)-ep(1,i))
411 . +lamda(i)*(depszz(i)-ep(3,i))
412 signzz(i)=sigozz(i)+wave(i)*(depszz(i)-ep(3,i))
413 . +lamda(i)*(depsxx(i)-ep(1,i))
414 . +lamda(i)*(depsyy(i)-ep(2,i))
415 signxy(i)=sigoxy(i)+g2(i)*(depsxy(i)-ep(4,i))
416 signyz(i)=sigoyz(i)+g2(i)*(depsyz(i)-ep(5,i))
417 signzx(i)=sigozx(i)+g2(i)*(depszx(i)-ep(6,i))
418c==========Module tangent=====================================
419
420 IF( dep(i) > zero)THEN
421 sigm(i)=-(signxx(i)+signyy(i)+signzz(i))*third
422 sigy=sqrt(three_half*(
423 . (signxx(i)+sigm(i))*(signxx(i)+sigm(i))
424 . +(signyy(i)+sigm(i))*(signyy(i)+sigm(i))
425 . +(signzz(i)+sigm(i))*(signzz(i)+sigm(i))
426 . +two*(signxy(i)*signxy(i)+signyz(i)*signyz(i)
427 . +signzx(i)*signzx(i))))
428 h = (sigy-uvar(i,4))/dep(i)
429 h = max(em10, h)
430 uvar(i,4) = sigy
431 uvar(i,5) = dep(i)
432 uvar(i,6) = max_asta(i)
433 etse(i) = h/g2(i)
434 pla(i) = pla(i)+dep(i)
435 ELSE
436 etse(i) = one
437 ENDIF
438C=======================================================================
439 soundsp(i) = sqrt((rbulk(i)+four_over_3*shear(i))/rho0(i))
440 yld(i) = uvar(i,4)
441 ENDDO
442c-----------
443 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
max
#define max(a, b)
Definition macros.h:21