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

Go to the source code of this file.

Functions/Subroutines

subroutine sigeps95 (nel, nuparam, nuvar, nfunc, ifunc, npf, tf, time, timestep, uparam, rho0, rho, volume, eint, ngl, 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, mfxx, mfxy, mfxz, mfyx, mfyy, mfyz, mfzx, mfzy, mfzz, soundsp, viscmax, uvar, off, ismstr, et, ihet, offg, epsth3, iexpan, nparf, uparamf, uvarf, nvarf, jcvt, gama_r)
subroutine rottoloc (nel, fp, r1x, r1y, r1z, r2x, r2y, r2z, r3x, r3y, r3z)
subroutine rottoglob (nel, fp, r1x, r1y, r1z, r2x, r2y, r2z, r3x, r3y, r3z)

Function/Subroutine Documentation

◆ rottoglob()

subroutine rottoglob ( integer, intent(in) nel,
dimension(nel,3,3), intent(inout) fp,
dimension(nel), intent(in) r1x,
dimension(nel), intent(in) r1y,
dimension(nel), intent(in) r1z,
dimension(nel), intent(in) r2x,
dimension(nel), intent(in) r2y,
dimension(nel), intent(in) r2z,
dimension(nel), intent(in) r3x,
dimension(nel), intent(in) r3y,
dimension(nel), intent(in) r3z )

Definition at line 492 of file sigeps95.F.

494C-----------------------------------------------
495C I m p l i c i t T y p e s
496C-----------------------------------------------
497#include "implicit_f.inc"
498C-----------------------------------------------
499C D u m m y A r g u m e n t s
500C-----------------------------------------------
501 INTEGER, INTENT(IN) :: NEL
502
503 my_real, DIMENSION(NEL), INTENT(IN) :: r1x(nel),r1y(nel),r1z(nel),
504 . r2x(nel), r2y(nel),r2z(nel),r3x(nel),r3y(nel),r3z(nel)
505
506 my_real, DIMENSION(NEL), INTENT(INOUT) :: fp(nel,3,3)
507C-----------------------------------------------
508C L o c a l V a r i a b l e s
509C-----------------------------------------------
510 INTEGER I
511 my_real
512 . sx,sy,sz,fxx,fxy,fxz,fyx,fyy,fyz,fzx,fzy,fzz
513C-----------------------------------------------
514 DO i=1,nel
515 sx = fp(i,1,1)*r1x(i)+fp(i,2,1)*r2x(i)+fp(i,3,1)*r3x(i)
516 sy = fp(i,1,2)*r1x(i)+fp(i,2,2)*r2x(i)+fp(i,3,2)*r3x(i)
517 sz = fp(i,1,3)*r1x(i)+fp(i,2,3)*r2x(i)+fp(i,3,3)*r3x(i)
518 fxx = sx*r1x(i)+sy*r2x(i)+sz*r3x(i)
519 fxy = sx*r1y(i)+sy*r2y(i)+sz*r3y(i)
520 fxz = sx*r1z(i)+sy*r2z(i)+sz*r3z(i)
521 sx = fp(i,1,1)*r1y(i)+fp(i,2,1)*r2y(i)+fp(i,3,1)*r3y(i)
522 sy = fp(i,1,2)*r1y(i)+fp(i,2,2)*r2y(i)+fp(i,3,2)*r3y(i)
523 sz = fp(i,1,3)*r1y(i)+fp(i,2,3)*r2y(i)+fp(i,3,3)*r3y(i)
524 fyx = sx*r1x(i)+sy*r2x(i)+sz*r3x(i)
525 fyy = sx*r1y(i)+sy*r2y(i)+sz*r3y(i)
526 fyz = sx*r1z(i)+sy*r2z(i)+sz*r3z(i)
527 sx = fp(i,1,1)*r1z(i)+fp(i,2,1)*r2z(i)+fp(i,3,1)*r3z(i)
528 sy = fp(i,1,2)*r1z(i)+fp(i,2,2)*r2z(i)+fp(i,3,2)*r3z(i)
529 sz = fp(i,1,3)*r1z(i)+fp(i,2,3)*r2z(i)+fp(i,3,3)*r3z(i)
530 fzx = sx*r1x(i)+sy*r2x(i)+sz*r3x(i)
531 fzy = sx*r1y(i)+sy*r2y(i)+sz*r3y(i)
532 fzz = sx*r1z(i)+sy*r2z(i)+sz*r3z(i)
533
534 fp(i,1,1)=fxx
535 fp(i,1,2)=fxy
536 fp(i,1,3)=fxz
537 fp(i,2,1)=fyx
538 fp(i,2,2)=fyy
539 fp(i,2,3)=fyz
540 fp(i,3,1)=fzx
541 fp(i,3,2)=fzy
542 fp(i,3,3)=fzz
543 ENDDO
544C-----------
545 RETURN
my_real
#define my_real
Definition cppsort.cpp:32

◆ rottoloc()

subroutine rottoloc ( integer, intent(in) nel,
dimension(nel,3,3), intent(inout) fp,
dimension(nel), intent(in) r1x,
dimension(nel), intent(in) r1y,
dimension(nel), intent(in) r1z,
dimension(nel), intent(in) r2x,
dimension(nel), intent(in) r2y,
dimension(nel), intent(in) r2z,
dimension(nel), intent(in) r3x,
dimension(nel), intent(in) r3y,
dimension(nel), intent(in) r3z )

Definition at line 432 of file sigeps95.F.

434C-----------------------------------------------
435C I m p l i c i t T y p e s
436C-----------------------------------------------
437#include "implicit_f.inc"
438C-----------------------------------------------
439C D u m m y A r g u m e n t s
440C-----------------------------------------------
441 INTEGER, INTENT(IN) :: NEL
442
443 my_real, DIMENSION(NEL), INTENT(IN) :: r1x(nel),r1y(nel),r1z(nel),
444 . r2x(nel), r2y(nel),r2z(nel),r3x(nel),r3y(nel),r3z(nel)
445
446 my_real, DIMENSION(NEL), INTENT(INOUT) :: fp(nel,3,3)
447C-----------------------------------------------
448C L o c a l V a r i a b l e s
449C-----------------------------------------------
450 INTEGER I
451 my_real
452 . sx,sy,sz,fxx,fxy,fxz,fyx,fyy,fyz,fzx,fzy,fzz
453C-----------------------------------------------
454 DO i=1,nel
455 sx = fp(i,1,1)*r1x(i)+fp(i,2,1)*r1y(i)+fp(i,3,1)*r1z(i)
456 sy = fp(i,1,2)*r1x(i)+fp(i,2,2)*r1y(i)+fp(i,3,2)*r1z(i)
457 sz = fp(i,1,3)*r1x(i)+fp(i,2,3)*r1y(i)+fp(i,3,3)*r1z(i)
458 fxx = sx*r1x(i)+sy*r1y(i)+sz*r1z(i)
459 fxy = sx*r2x(i)+sy*r2y(i)+sz*r2z(i)
460 fxz = sx*r3x(i)+sy*r3y(i)+sz*r3z(i)
461 sx = fp(i,1,1)*r2x(i)+fp(i,2,1)*r2y(i)+fp(i,3,1)*r2z(i)
462 sy = fp(i,1,2)*r2x(i)+fp(i,2,2)*r2y(i)+fp(i,3,2)*r2z(i)
463 sz = fp(i,1,3)*r2x(i)+fp(i,2,3)*r2y(i)+fp(i,3,3)*r2z(i)
464 fyx = sx*r1x(i)+sy*r1y(i)+sz*r1z(i)
465 fyy = sx*r2x(i)+sy*r2y(i)+sz*r2z(i)
466 fyz = sx*r3x(i)+sy*r3y(i)+sz*r3z(i)
467 sx = fp(i,1,1)*r3x(i)+fp(i,2,1)*r3y(i)+fp(i,3,1)*r3z(i)
468 sy = fp(i,1,2)*r3x(i)+fp(i,2,2)*r3y(i)+fp(i,3,2)*r3z(i)
469 sz = fp(i,1,3)*r3x(i)+fp(i,2,3)*r3y(i)+fp(i,3,3)*r3z(i)
470 fzx = sx*r1x(i)+sy*r1y(i)+sz*r1z(i)
471 fzy = sx*r2x(i)+sy*r2y(i)+sz*r2z(i)
472 fzz = sx*r3x(i)+sy*r3y(i)+sz*r3z(i)
473 fp(i,1,1)=fxx
474 fp(i,1,2)=fxy
475 fp(i,1,3)=fxz
476 fp(i,2,1)=fyx
477 fp(i,2,2)=fyy
478 fp(i,2,3)=fyz
479 fp(i,3,1)=fzx
480 fp(i,3,2)=fzy
481 fp(i,3,3)=fzz
482 ENDDO
483C-----------
484 RETURN

◆ sigeps95()

subroutine sigeps95 ( integer, intent(in) nel,
integer, intent(in) nuparam,
integer, intent(in) nuvar,
integer nfunc,
integer, dimension(nfunc) ifunc,
integer, dimension(*) npf,
tf,
time,
timestep,
intent(in) uparam,
dimension(nel) rho0,
dimension(nel) rho,
dimension(nel) volume,
dimension(nel) eint,
integer, dimension(nel), intent(in) ngl,
dimension(nel) epspxx,
dimension(nel) epspyy,
dimension(nel) epspzz,
dimension(nel) epspxy,
dimension(nel) epspyz,
dimension(nel) epspzx,
dimension(nel) depsxx,
dimension(nel) depsyy,
dimension(nel) depszz,
dimension(nel) depsxy,
dimension(nel) depsyz,
dimension(nel) depszx,
dimension(nel) epsxx,
dimension(nel) epsyy,
dimension(nel) epszz,
dimension(nel) epsxy,
dimension(nel) epsyz,
dimension(nel) epszx,
dimension(nel) sigoxx,
dimension(nel) sigoyy,
dimension(nel) sigozz,
dimension(nel) sigoxy,
dimension(nel) sigoyz,
dimension(nel) sigozx,
signxx,
signyy,
signzz,
signxy,
signyz,
signzx,
sigvxx,
sigvyy,
sigvzz,
sigvxy,
sigvyz,
sigvzx,
dimension(nel) mfxx,
dimension(nel) mfxy,
dimension(nel) mfxz,
dimension(nel) mfyx,
dimension(nel) mfyy,
dimension(nel) mfyz,
dimension(nel) mfzx,
dimension(nel) mfzy,
dimension(nel) mfzz,
soundsp,
viscmax,
uvar,
off,
integer, intent(in) ismstr,
et,
integer, intent(in) ihet,
dimension(nel) offg,
dimension(nel) epsth3,
integer, intent(in) iexpan,
integer, intent(in) nparf,
intent(in) uparamf,
uvarf,
integer, intent(in) nvarf,
integer, intent(in) jcvt,
gama_r )

Definition at line 40 of file sigeps95.F.

55C-----------------------------------------------
56C I M P L I C I T T Y P E S
57C-----------------------------------------------
58#include "implicit_f.inc"
59C-----------------------------------------------
60C C O M M O N
61C-----------------------------------------------
62#include "impl1_c.inc"
63C----------------------------------------------------------------
64C I N P U T A R G U M E N T S
65C----------------------------------------------------------------
66 INTEGER ,INTENT(IN) :: NPARF
67 INTEGER ,INTENT(IN) :: NEL,JCVT,NVARF,NUPARAM,NUVAR,ISMSTR,IHET,IEXPAN
68 INTEGER ,INTENT(IN) :: NGL(NEL)
69 my_real :: time,timestep
70 my_real ,DIMENSION(NUPARAM) ,INTENT(IN) :: uparam
71 my_real ,DIMENSION(NPARF) ,INTENT(IN) :: uparamf
72 my_real ,DIMENSION(NEL) :: rho,rho0,volume,eint,epsth3,
73 . epspxx, epspyy, epspzz, epspxy, epspyz, epspzx,
74 . depsxx, depsyy, depszz, depsxy, depsyz, depszx,
75 . epsxx , epsyy , epszz , epsxy , epsyz , epszx ,
76 . sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx,offg,
77 . mfxx,mfxy,mfxz,mfyx,mfyy,mfyz,mfzx,mfzy,mfzz
78
79 my_real :: gama_r(nel,6)
80C----------------------------------------------------------------
81C O U T P U T A R G U M E N T S
82C----------------------------------------------------------------
83 my_real
84 . signxx(nel), signyy(nel), signzz(nel),
85 . signxy(nel), signyz(nel), signzx(nel),
86 . sigvxx(nel), sigvyy(nel), sigvzz(nel),
87 . sigvxy(nel), sigvyz(nel), sigvzx(nel),
88 . soundsp(nel), viscmax(nel), et(nel)
89C----------------------------------------------------------------
90C I N P U T O U T P U T A R G U M E N T S
91C----------------------------------------------------------------
92 my_real
93 . uvar(nel,nuvar), off(nel) ,uvarf(nel,nvarf)
94C----------------------------------------------------------------
95C VARIABLES FOR FUNCTION INTERPOLATION
96C----------------------------------------------------------------
97 INTEGER NPF(*), NFUNC, IFUNC(NFUNC)
98 my_real finter,fintte,tf(*),fint2v
99 EXTERNAL finter,fintte
100C----------------------------------------------------------------
101C L O C A L V A R I B L E S
102C----------------------------------------------------------------
103 INTEGER I,J,K,KK,LL,FLAGBB,DIRECT,ITER,NITER,FLAG_MUL,IAVIS,IFORM
104
105 my_real
106 . et1,et2,et3,a1,a10,ksi,mu,d,lam,g,rbulk,aa,bb,cc,sb, temp,
107 . temp2,maxl,stiff0,dsig,deps,coef1,coef2,coef3,coef4,coef5,coef6,
108 . aa1,b1,c1,dd1,aa2,b2,c2,dd2,aa3,b3,c3,dd3,inv12,inv22,inv32,
109 . evb12,evb22,evb32,jdet2,jm1,etvol,factor,tauref,r3r3,r2r2,r1r1,
110 . c10,c01,c20,c11,c02,c30,c21,c12,c03,d1,d2,d3,expc,expm,r(3),jcst(3),
111 . bi1(nel),bi2(nel),jdet(nel),stiff(nel),
112 . ta1(nel), ta2(nel),ta3(nel),t1(nel), t2(nel),t3(nel),lpchain(nel),
113 . tb1(nel), tb2(nel),tb3(nel),trace(nel),traceb(nel),eta(nel),etb(nel),
114 . sb1(nel), sb2(nel),sb3(nel),sb4(nel), sb5(nel),sb6(nel),
115 . tbnorm(nel),dgamma(nel),
116 . devb1(nel),devb2(nel),devb3(nel),devb4(nel),devb5(nel),devb6(nel),
117 . r1x(nel),r1y(nel),r1z(nel),r2x(nel),r2y(nel),r2z(nel),r3x(nel),r3y(nel),r3z(nel),
118 . nb(nel,3),
119 . f(nel,3,3),ft(nel,3,3),fe(nel,3,3),invfe(nel,3,3),fet(nel,3,3),fp(nel,3,3),
120 . fft(nel,3,3),invfpo(nel,3,3),matb(nel,3,3),fpo(nel,3,3),bpo(nel,3,3),invsn(nel,3,3),
121 . sig(nel,3,3),sigb(nel,3,3),siga(nel,3,3),sn(nel,3,3),fedp(nel,3,3),
122 . dfp(nel,3,3),lb(nel,3,3),dfp2(nel,3,3),fpdot(nel,3,3),
123 . caii(nel,3),cbii(nel,3),c_max(nel),cii(3)
124 my_real
125 . coefr,betaf ,coefm
126C----------------------------------------------------------------
127 coefr = one
128 betaf = zero
129 coefm = one
130
131c material parameters
132 c10 = uparam(1)
133 c01 = uparam(2)
134 c20 = uparam(3)
135 c11 = uparam(4)
136 c02 = uparam(5)
137 c30 = uparam(6)
138 c21 = uparam(7)
139 c12 = uparam(8)
140 c03 = uparam(9)
141 d1 = uparam(11) !ONE/D1
142 d2 = uparam(12) !ONE/D2
143 d3 = uparam(13) !ONE/D3
144 sb = uparam(14)
145 a10 = uparam(15)
146 a1 = a10*timestep
147 expc = uparam(16)
148 expm = uparam(17)
149 ksi = uparam(18)
150 g = uparam(19)
151 rbulk= uparam(20)
152 tauref= uparam(22)
153 iform = nint(uparam(23))
154 !
155 stiff0 = four_over_3*g + rbulk
156 maxl = one
157 iavis =1
158 IF (a10*sb==zero) iavis=0
159
160 flag_mul = uparam(21) !CALCULATED IN UPDMAT = 1 IF IRUP==33
161 IF (flag_mul > zero)THEN
162 coefr = uparamf(1)
163 betaf = uparamf(2)
164 coefm = uparamf(3)
165 ENDIF
166
167C
168 IF (iavis>0) THEN
169 DO i = 1, nel
170 !retrieve previous viscous deformation gradient
171 fpo(i,1,1) = uvar(i,1)
172 fpo(i,2,2) = uvar(i,2)
173 fpo(i,3,3) = uvar(i,3)
174 fpo(i,1,2) = uvar(i,4)
175 fpo(i,2,3) = uvar(i,5)
176 fpo(i,3,1) = uvar(i,6)
177 fpo(i,2,1) = uvar(i,7)
178 fpo(i,3,2) = uvar(i,8)
179 fpo(i,1,3) = uvar(i,9)
180 ENDDO
181
182 IF (jcvt >0 ) THEN ! corotational => need to rotate Fp in the LOCAL frame
183 r1x(1:nel) = gama_r(1:nel,1)
184 r1y(1:nel) = gama_r(1:nel,2)
185 r1z(1:nel) = gama_r(1:nel,3)
186 r2x(1:nel) = gama_r(1:nel,4)
187 r2y(1:nel) = gama_r(1:nel,5)
188 r2z(1:nel) = gama_r(1:nel,6)
189
190 r3x(1:nel) = r1y(1:nel)* r2z(1:nel) - r1z(1:nel) * r2y(1:nel)
191 r3y(1:nel) = r1z(1:nel)* r2x(1:nel) - r1x(1:nel) * r2z(1:nel)
192 r3z(1:nel) = r1x(1:nel)* r2y(1:nel) - r1y(1:nel) * r2x(1:nel)
193 DO i=1,nel
194 r3r3 = sqrt(r3x(i)*r3x(i) + r3y(i)*r3y(i) + r3z(i)*r3z(i))
195 IF (r3r3 /= zero) THEN
196 r3x(i) = r3x(i)/r3r3
197 r3y(i) = r3y(i)/r3r3
198 r3z(i) = r3z(i)/r3r3
199 ENDIF
200 ENDDO
201 CALL rottoloc (nel,fpo,
202 . r1x, r1y, r1z, r2x, r2y, r2z, r3x, r3y, r3z)
203
204 ENDIF
205 ENDIF !(IAVIS>0)
206
207 DO i=1,nel
208 f(i,1,1) = one+mfxx(i)
209 f(i,2,2) = one+mfyy(i)
210 f(i,3,3) = one+mfzz(i)
211 f(i,1,2) = mfxy(i)
212 f(i,2,3) = mfyz(i)
213 f(i,3,1) = mfzx(i)
214 f(i,2,1) = mfyx(i)
215 f(i,3,2) = mfzy(i)
216 f(i,1,3) = mfxz(i)
217 ENDDO
218
219
220 CALL prodaat(f , fft, nel) ! b = F * FT
221
222
223 !-------------------------
224 !chain A: compute stress : cauchy stress in chain A
225 !-------------------------
226 !entry B=F*FT matrix -- output hyperelastic stress
227 CALL polystrest2(
228 1 nel , fft , c10 , c01, c20,
229 2 c11 ,c02 , c30 , c21, c12,
230 3 c03 ,d1 , d2 , d3, siga ,
231 4 bi1 ,bi2 ,jdet ,flag_mul,
232 5 nvarf,coefr,betaf,coefm ,
233 6 uvarf,caii ,rbulk,iform)
234
235
236 IF (iavis>0) THEN
237 !-------------------------
238 !chain B: compute stress :
239 !-------------------------
240 DO i=1,nel
241 sigb(i,1,1) = zero
242 sigb(i,2,2) = zero
243 sigb(i,3,3) = zero
244 sigb(i,1,2) = zero
245 sigb(i,2,3) = zero
246 sigb(i,3,1) = zero
247 ENDDO
248
249 !CALL CALCMATB (NEL, F, FPO, MATB)
250 !FE ={F}{Fp_old}^(-1) then MATB = FE FE^(T) (elastic part)
251 CALL kmatinv3 (fpo , invfpo, nel) !INVFPO = INVERSE (FP)
252 CALL prodmat (f , invfpo, fe, nel) ! FE = F * INVFPO
253 CALL prodaat (fe , matb , nel) ! MATB = FE FE^(T)
254
255 !---------------------------------
256 !ESTIMATE TRIAL CAUCHY STRESS : cauchy stress in chain B
257 !---------------------------------
258 !entry B matrix -- output hyperelastic stress
259 CALL polystress2(
260 1 nel , matb , c10, c01, c20,
261 2 c11 ,c02, c30, c21, c12,
262 3 c03 ,d1 ,d2 , d3, sigb ,
263 4 bi1,bi2,jdet ,flag_mul,
264 5 nvarf,coefr, betaf,coefm ,
265 6 uvarf,rbulk,iform)
266
267
268 DO i=1,nel ! scale trial cauchy stress in chain B
269 sigb(i,1,1) = sb * sigb(i,1,1)
270 sigb(i,2,2) = sb * sigb(i,2,2)
271 sigb(i,3,3) = sb * sigb(i,3,3)
272 sigb(i,1,2) = sb * sigb(i,1,2)
273 sigb(i,2,3) = sb * sigb(i,2,3)
274 sigb(i,3,1) = sb * sigb(i,3,1)
275 ENDDO
276
277
278 DO i=1,nel
279 traceb(i) = third*(sigb(i,1,1) +sigb(i,2,2) + sigb(i,3,3))
280 !DEVIATOR OF STRESS B
281 sb1(i) = sigb(i,1,1) - traceb(i)
282 sb2(i) = sigb(i,2,2) - traceb(i)
283 sb3(i) = sigb(i,3,3) - traceb(i)
284 sb4(i) = sigb(i,1,2)
285 sb5(i) = sigb(i,2,3)
286 sb6(i) = sigb(i,3,1)
287 !Nomr of stress B
288 tbnorm(i) = sqrt(max(em20,sb1(i)**2+sb2(i)**2+sb3(i)**2
289 . + two*(sb4(i)**2+sb5(i)**2+sb6(i)**2 )) ) ! NORM!
290 ENDDO
291
292 !COMPUTE EEFFECTIVE CREEP STRAIN RATE
293 !------------------------------------
294 !DGAMMA(I) = A1* ((LPCHAIN(I) - ONE +KSI)**C) * (TBNORM(I)/TAUREF)**EXPM
295 CALL viscbb ( nel , fpo, tbnorm , a1 , expc,
296 . expm, ksi, tauref , dgamma )
297
298
299 IF (time == zero) dgamma(1:nel) = uvar(1:nel,10)
300 DO i=1,nel ! calcul lde Lb
301 uvar(i,10)= dgamma(i)
302 factor = dgamma(i)/tbnorm(i)
303 lb(i,1,1) = factor *sb1(i) !gamma_dot*N
304 lb(i,2,2) = factor *sb2(i)
305 lb(i,3,3) = factor *sb3(i)
306 lb(i,1,2) = factor *sb4(i) !gamma_dot*N
307 lb(i,2,3) = factor *sb5(i)
308 lb(i,3,1) = factor *sb6(i)
309 lb(i,2,1) = lb(i,1,2)
310 lb(i,3,2) = lb(i,2,3)
311 lb(i,1,3) = lb(i,3,1)
312
313 ENDDO
314 !------------------------------------
315 !Solve F_n+1 viscous :
316 !------------------------------------
317 !PRODMAT(A, B, C, NEL) computes C(NEL,3,3) which is the product [C] = [A][B]
318 CALL kmatinv3(fe , invfe, nel) ! INVFE = INVERSE (FE)
319 CALL prodmat (lb , fe, fedp, nel) ! FEDP = Lb * FE
320 CALL prodmat(invfe ,fedp, dfp , nel) ! DFP = INVFE * Lb * FE
321 !not needed CALL PRODMAT(DFP ,FPO, FPDOT , NEL) !FPDOT = DFP * FPO
322 !CALL PRODMAT(DFP ,DFP, DFP2 , NEL) ! if order 2
323 DO i=1,nel
324 ! Fp_n+1 = exp(DFP) * Fp_old = (I+DFP) * Fp_old
325 sn(i,1,1) = one + dfp(i,1,1) !+ HALF * DFP2(I,1,1)
326 sn(i,2,2) = one + dfp(i,2,2) !+ HALF * DFP2(I,2,2)
327 sn(i,3,3) = one + dfp(i,3,3) !+ HALF * DFP2(I,3,3)
328 sn(i,1,2) = dfp(i,1,2)!+ HALF * DFP2(I,1,2)
329 sn(i,2,3) = dfp(i,2,3)!+ HALF * DFP2(I,2,3)
330 sn(i,3,1) = dfp(i,3,1)!+ HALF * DFP2(I,3,1)
331 sn(i,2,1) = dfp(i,2,1)!+ HALF * DFP2(I,2,1)
332 sn(i,3,2) = dfp(i,3,2)!+ HALF * DFP2(I,3,2)
333 sn(i,1,3) = dfp(i,1,3)!+ HALF * DFP2(I,1,3)
334 ENDDO
335 CALL prodmat(sn ,fpo, fp, nel) ! Fp_N+1 = (I + dT *DFP)* Fp_N
336 CALL calcmatb (nel, f, fp, matb) !FE ={F}{Fp}^(-1) then MATB = FE FE^(T)
337
338 !-------------------------
339 !chain B: compute stress : cauchy stress in chain B
340 !-------------------------
341 CALL polystrest2(
342 1 nel , matb , c10, c01, c20,
343 2 c11 ,c02, c30, c21, c12,
344 3 c03 ,d1 ,d2 , d3, sigb ,
345 4 bi1,bi2,jdet ,flag_mul,
346 5 nvarf,coefr, betaf,coefm ,
347 6 uvarf,cbii,rbulk,iform)
348
349 DO i=1,nel ! scale cauchy stress in chain B
350 sigb(i,1,1) = sb * sigb(i,1,1)
351 sigb(i,2,2) = sb * sigb(i,2,2)
352 sigb(i,3,3) = sb * sigb(i,3,3)
353 sigb(i,1,2) = sb * sigb(i,1,2)
354 sigb(i,2,3) = sb * sigb(i,2,3)
355 sigb(i,3,1) = sb * sigb(i,3,1)
356 ENDDO
357 IF (jcvt >0 ) THEN ! corotational => need to rotate NEW Fp TO the GLOBAL frame
358 CALL rottoglob (nel,fp,
359 . r1x, r1y, r1z, r2x, r2y, r2z, r3x, r3y, r3z)
360
361 ENDIF
362 DO i=1,nel
363 uvar(i,1) = fp(i,1,1)!store visous gradient of deformation
364 uvar(i,2) = fp(i,2,2)
365 uvar(i,3) = fp(i,3,3)
366 uvar(i,4) = fp(i,1,2)
367 uvar(i,5) = fp(i,2,3)
368 uvar(i,6) = fp(i,3,1)
369 uvar(i,7) = fp(i,2,1)
370 uvar(i,8) = fp(i,3,2)
371 uvar(i,9) = fp(i,1,3)
372 ENDDO
373C----- A10*SB=0 ->sigB=sigA
374 ELSE
375 DO i=1,nel ! scale cauchy stress in chain B
376 sigb(i,1:3,1:3) = sb * siga(i,1:3,1:3)
377 cbii(i,1:3) = caii(i,1:3)
378 ENDDO
379 END IF !(IAVIS>0) THEN
380
381 !--------------------
382 !UPDATE TOTAL STRESS : sum of the stresses of each network
383 !--------------------
384 DO i=1,nel
385 signxx(i) = siga(i,1,1) + sigb(i,1,1)
386 signyy(i) = siga(i,2,2) + sigb(i,2,2)
387 signzz(i) = siga(i,3,3) + sigb(i,3,3)
388 signxy(i) = siga(i,1,2) + sigb(i,1,2)
389 signyz(i) = siga(i,2,3) + sigb(i,2,3)
390 signzx(i) = siga(i,3,1) + sigb(i,3,1)
391 ENDDO
392
393
394 DO i=1,nel
395 cii(1:3)=caii(i,1:3)+sb*cbii(i,1:3)
396 c_max(i)=max(stiff0,cii(1),cii(2),cii(3))
397 ENDDO
398
399 DO i=1,nel
400c DSIG = SQRT((SIGNXX(I)-SIGOXX(I))**2+(SIGNYY(I)-SIGOYY(I))**2
401c . +(SIGNZZ(I)-SIGOZZ(I))**2
402c . + TWO*((SIGNXY(I)-SIGOXY(I))**2
403c . +(SIGNYZ(I)-SIGOYZ(I))**2+(SIGNZX(I)-SIGOZX(I))**2))
404c
405c DEPS = SQRT (DEPSXX(I)**2 + DEPSYY(I)**2 + DEPSZZ(I)**2
406c . + TWO*(DEPSXY(I)**2 + DEPSYZ(I)**2 + DEPSZX(I)**2))
407c
408c IF(DEPS == ZERO)THEN
409c STIFF(I)= STIFF0
410c ELSE
411c STIFF(I)= MAX(STIFF0,DSIG /MAX(EM20,DEPS))
412c ENDIF
413C
414 soundsp(i)=sqrt(c_max(i)/rho(i))
415 viscmax(i) = zero
416 ENDDO
417 IF (impl_s > 0 .OR. ihet > 1) THEN
418 DO i=1,nel
419 et(i)= one
420 ENDDO
421 ENDIF
422
423C
424 RETURN
calcmatb
subroutine calcmatb(nel, f, fp, matb)
Definition calcmatb.F:39
kmatinv3
subroutine kmatinv3(mat, ainv, nel)
Definition kmatinv.F:32
max
#define max(a, b)
Definition macros.h:21
prodaat
subroutine prodaat(a, c, nel)
Definition prodAAT.F:34
prodmat
subroutine prodmat(a, b, c, nel)
Definition prodmat.F:35
rottoloc
subroutine rottoloc(nel, fp, r1x, r1y, r1z, r2x, r2y, r2z, r3x, r3y, r3z)
Definition sigeps95.F:434
rottoglob
subroutine rottoglob(nel, fp, r1x, r1y, r1z, r2x, r2y, r2z, r3x, r3y, r3z)
Definition sigeps95.F:494
polystress2
subroutine polystress2(nel, matb, c10, c01, c20, c11, c02, c30, c21, c12, c03, d1, d2, d3, sig, bi1, bi2, jdet, flag_mul, nvarf, coefr, betaf, coefm, uvarf, rbulk, iform)
Definition sigpoly.F:43
polystrest2
subroutine polystrest2(nel, matb, c10, c01, c20, c11, c02, c30, c21, c12, c03, d1, d2, d3, sig, bi1, bi2, jdet, flag_mul, nvarf, coefr, betaf, coefm, uvarf, cii, rbulk, iform)
Definition sigpoly.F:204
viscbb
subroutine viscbb(nel, fp, tbnorm, a1, expc, expm, ksi, tauref, dgamma)
Definition viscbb.F:33