OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps42.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!|| sigeps42 ../engine/source/materials/mat/mat042/sigeps42.F
25!||--- called by ------------------------------------------------------
26!|| mulaw ../engine/source/materials/mat_share/mulaw.F90
27!|| mulaw8 ../engine/source/materials/mat_share/mulaw8.F90
28!||--- calls -----------------------------------------------------
29!|| finter ../engine/source/tools/curve/finter.F
30!|| prodaat ../engine/source/materials/tools/prodAAT.F
31!|| prodmat ../engine/source/materials/tools/prodmat.F
32!|| valpvec_v ../engine/source/materials/mat/mat033/sigeps33.F
33!|| valpvecdp_v ../engine/source/materials/mat/mat033/sigeps33.F
34!||====================================================================
35 SUBROUTINE sigeps42(
36 1 NEL ,NUPARAM ,NUVAR ,NFUNC ,IFUNC ,NPF ,
37 2 TF ,TIME ,TIMESTEP,UPARAM ,RHO0 ,RHO ,
38 3 VOLUME ,EINT ,UVAR ,OFF ,OFFG ,SOUNDSP ,
39 4 EPSP1 ,EPSP2 ,EPSP3 ,EPSP4 ,EPSP5 ,EPSP6 ,
40 5 EPSXX ,EPSYY ,EPSZZ ,EPSXY ,EPSYZ ,EPSZX ,
41 6 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
42 7 MFXX ,MFXY ,MFXZ ,MFYX ,MFYY ,MFYZ ,
43 8 MFZX ,MFZY ,MFZZ ,VISCMAX ,ISMSTR ,ET ,
44 9 IHET ,EPSTH3 ,IEXPAN ,NIPARAM ,IPARAM )
45C-----------------------------------------------
46C I M P L I C I T T Y P E S
47C-----------------------------------------------
48#include "implicit_f.inc"
49C-----------------------------------------------
50C G L O B A L P A R A M E T E R S
51C-----------------------------------------------
52#include "mvsiz_p.inc"
53C-----------------------------------------------
54C C O M M O N
55C-----------------------------------------------
56#include "scr05_c.inc"
57#include "scr17_c.inc"
58#include "impl1_c.inc"
59#include "tabsiz_c.inc"
60C----------------------------------------------------------------
61C I N P U T A R G U M E N T S
62C----------------------------------------------------------------
63 INTEGER, INTENT(IN) :: NUPARAM
64 INTEGER, INTENT(IN) :: NIPARAM
65 INTEGER, INTENT(IN) :: NEL,NUVAR,ISMSTR,IHET,IEXPAN
66 my_real, INTENT(IN) :: TIME, TIMESTEP
67 INTEGER, DIMENSION(NIPARAM), INTENT(IN) :: IPARAM(NIPARAM)
68 my_real, DIMENSION(NUPARAM), INTENT(IN) :: UPARAM(NUPARAM)
69 my_real, DIMENSION(NEL), INTENT(IN) ::RHO,RHO0,VOLUME,EINT,
70 . OFFG,EPSTH3,EPSP1,EPSP2,EPSP3,EPSP4,EPSP5,EPSP6,
71 . EPSXX,EPSYY,EPSZZ,EPSXY,EPSYZ,EPSZX,
72 . mfxx,mfxy,mfxz,mfyx,mfyy,mfyz,mfzx,mfzy,mfzz
73C----------------------------------------------------------------
74C O U T P U T A R G U M E N T S
75C----------------------------------------------------------------
76 my_real ,DIMENSION(NEL), INTENT(OUT) :: soundsp,viscmax,
77 . signxx,signyy,signzz,signxy,signyz,signzx
78C----------------------------------------------------------------
79C I N P U T O U T P U T A R G U M E N T S
80C----------------------------------------------------------------
81 my_real, DIMENSION(NEL,NUVAR), INTENT(INOUT) :: uvar(nel,nuvar)
82 my_real, DIMENSION(NEL), INTENT(INOUT) :: off(nel),et(nel)
83C----------------------------------------------------------------
84C VARIABLES FOR FUNCTION INTERPOLATION
85C----------------------------------------------------------------
86 INTEGER NPF(SNPC),NFUNC,IFUNC(NFUNC)
87 my_real FINTER,TF(STF)
88 EXTERNAL FINTER
89C----------------------------------------------------------------
90C L O C A L V A R I A B L E S
91C----------------------------------------------------------------
92 INTEGER I,II,J,K,KFP,NPRONY,NORDER,IFORM
93 my_real TENSIONCUT,AMAX,GVMAX,GMAX,FFAC,TRACE,RBULK,DPDMU,
94 . AA,CC,FAC,INVE,ETV,RV_PUI,AMIN,NU_1
95 my_real, DIMENSION(3) :: lam_al,fft
96 my_real, DIMENSION(10) :: al,mu
97 my_real, DIMENSION(100) :: gi,taux
98 my_real, DIMENSION(NEL) :: sumdwdl,rv,j2third,p,dwdrv,eti,gtmax
99 my_real, DIMENSION(3,3) :: a
100 my_real, DIMENSION(NEL,3) :: t1,sigprv,ev,evm,dwdl,cii
101 my_real, DIMENSION(NEL,6) :: dotb,svisc,c0,c1
102 my_real, DIMENSION(NEL,3,3) :: f,ll,bb,lb,blt
103 my_real, DIMENSION(NEL,6,100) :: h0, h
104 my_real, DIMENSION(MVSIZ,3) :: evv
105 my_real, DIMENSION(MVSIZ,6) :: av
106 my_real, DIMENSION(MVSIZ,3,3) :: dirprv
107 my_real, DIMENSION(NEL) :: p_fac
108 DOUBLE PRECISION AMAX1
109 !=======================================================================
110c
111 !=======================================================================
112 ! - INITIALISATION OF COMPUTATION ON TIME STEP
113 !=======================================================================
114 ! Recovering model parameters
115!
116 norder = iparam(1) ! Order of Ogden model
117 nprony = iparam(2) ! Number of Prony Series
118 iform = iparam(3) ! Flag for strain energy density formulation
119
120 gmax = zero
121 DO i=1,norder
122 ! -> Shear hyperelastic modulus parameters
123 mu(i) = uparam(i)
124 ! -> Material exponents
125 al(i) = uparam(10+i)
126 ! -> Shear modulus
127 gmax = gmax + mu(i)*al(i)
128 ENDDO
129 ! Bulk modulus
130 rbulk = uparam(21)
131 ! Cutoff stress in tension
132 tensioncut = uparam(23)
133 ! Bulk function scale factor
134 ffac = uparam(24)
135 ! Tabulated bulk function ID
136 kfp = ifunc(1)
137 ! Prony series parameters
138 gvmax = zero
139 etv = zero
140 IF (nprony > 0) THEN
141 DO i=1,nprony
142 taux(i) = uparam(24 + nprony + i)
143 gi(i) = uparam(24 + i)
144 gvmax = gvmax + gi(i)
145 ENDDO
146 etv = min(gvmax,rbulk)/gmax
147 ENDIF
148c
149 ! Fill strain tensor
150 DO i=1,nel
151 av(i,1) = epsxx(i)
152 av(i,2) = epsyy(i)
153 av(i,3) = epszz(i)
154 av(i,4) = epsxy(i) * half
155 av(i,5) = epsyz(i) * half
156 av(i,6) = epszx(i) * half
157 ENDDO
158c
159 ! Compute principal strains
160 ! Note: in simple precision, principal strains are computed
161 ! with double precision
162 IF (iresp == 1) THEN
163 CALL valpvecdp_v(av,evv,dirprv,nel)
164 ELSE
165 CALL valpvec_v(av,evv,dirprv,nel)
166 ENDIF
167c
168 ! Compute principal stretches depending on strain formulation
169 ! (Principal stretches = lambda_i)
170 ! -> Logarithmic strains
171 IF (ismstr == 0 .OR. ismstr == 2 .OR. ismstr == 4) THEN
172 DO i=1,nel
173 ev(i,1) = exp(evv(i,1))
174 ev(i,2) = exp(evv(i,2))
175 ev(i,3) = exp(evv(i,3))
176 ENDDO
177 ! -> Green-Lagrange strains
178 ELSEIF (ismstr == 10 .OR. ismstr == 12) THEN
179 DO i =1,nel
180 IF (offg(i)<=one) THEN
181 ev(i,1) = sqrt(evv(i,1) + one)
182 ev(i,2) = sqrt(evv(i,2) + one)
183 ev(i,3) = sqrt(evv(i,3) + one)
184 ELSE
185 ev(i,1) = evv(i,1) + one
186 ev(i,2) = evv(i,2) + one
187 ev(i,3) = evv(i,3) + one
188 END IF
189 ENDDO
190 ! -> Engineering strains
191 ELSE
192 DO i=1,nel
193 ev(i,1) = evv(i,1) + one
194 ev(i,2) = evv(i,2) + one
195 ev(i,3) = evv(i,3) + one
196 ENDDO
197 ENDIF
198c
199 ! Implicit simulation - Hourglass for HEPH
200 IF (impl_s > 0 .OR. ihet > 1) THEN
201 DO j = 1,3
202 ! ETI = sum[MU(i)*AL(i)*AMAX**(AL(i)-1)] (i=1,norder)
203 ! AMAX = 0 --> ETI = 0
204 ! else --> ETI = sum[MU(i)*AL(i)*exp( (AL(i)-1)*ln(AMAX) ) ] ; (i=1,norder)
205 eti(1:nel) = zero
206 DO ii = 1,norder
207 IF(mu(ii)*al(ii)/=zero) THEN
208 DO i=1,nel
209 amax = ev(i,j)
210 IF(amax/=zero) THEN
211 IF((al(ii)-one)/=zero) THEN
212 eti(i) = eti(i) + mu(ii)*al(ii) * exp((al(ii)-one)*log(amax))
213 ELSE
214 eti(i) = eti(i) + mu(ii)*al(ii)
215 ENDIF
216 ENDIF
217 ENDDO
218 ENDIF
219 ENDDO
220 et(1:nel)= max(eti(1:nel),et(1:nel))
221 ENDDO
222 DO i=1,nel
223 et(i) = max(one,et(i)/gmax)
224 et(i)= et(i)+etv
225 ENDDO
226 ENDIF
227c
228 ! Relative volume computation
229 ! (RHO/RHO0) = def(F) with F = Grad(Strain)
230 DO i=1,nel
231 rv(i) = (ev(i,1)*ev(i,2)*ev(i,3))
232 ENDDO
233
234 ! Considering thermal expansion
235 IF (iexpan > 0 .AND. (ismstr == 10 .OR. ismstr == 11 .OR. ismstr == 12)) THEN
236 DO i=1,nel
237 rv(i) = rv(i) - epsth3(i)
238 ENDDO
239 ENDIF
240c--- avoid buckling
241 p_fac(1:nel) = one
242 IF (rbulk > 24*gmax) THEN
243 nu_1 = fourty*(half-(3*rbulk-gmax)/(6*rbulk+gmax))
244 DO i=1,nel
245 amin = min(ev(i,1),ev(i,2),ev(i,3))
246 IF (amin<zep2) p_fac(i) = max(one,nu_1/max(em20,amin))
247 ENDDO
248 END IF
249c
250 ! Compute normalized (deviatoric) stretches and bulk modulus
251 ! (Deviatoric principal stretches = lambda_bar_i)
252 DO i=1,nel
253c
254 ! Deviatoric stretches computed only if relative volume is positive
255 IF (rv(i) > zero) THEN
256 rv_pui = exp((-third)*log(rv(i))) ! -> J^(-1/3)
257 j2third(i) = rv_pui**2 ! -> (J^(-1/3))^2
258 ELSE
259 rv_pui = zero
260 j2third(i) = zero
261 ENDIF
262 evm(i,1) = ev(i,1)*rv_pui
263 evm(i,2) = ev(i,2)*rv_pui
264 evm(i,3) = ev(i,3)*rv_pui
265c
266 ! Tabulated bulk modulus with respect to relative volume
267 IF (kfp > 0) THEN
268 p(i) = rbulk*ffac*finter(kfp,rv(i),npf,tf,dpdmu)
269 ! Constant bulk modulus
270 ELSE
271 p(i) = p_fac(i)*rbulk
272 ENDIF
273 ENDDO
274c
275 ! For HEPH Hourglass treatment and soundspeed computation
276 gtmax(1:nel) = gmax + gvmax
277 eti(1:nel) = zero
278 cii(1:nel,1:3) = zero
279 DO ii = 1,norder
280 IF (mu(ii)*al(ii) /= zero) THEN
281 DO i=1,nel
282 lam_al(1:3) = exp(al(ii)*log(evm(i,1:3)))
283 amax = third*(lam_al(1)+lam_al(2)+lam_al(3))
284 cii(i,1:3) = cii(i,1:3) + mu(ii)*al(ii)*(lam_al(1:3)+amax)
285 ENDDO
286 ENDIF
287 ENDDO
288 ! Note: 2GT=3/2 Cii (factor 3 already inside CII)
289 amax1 = 0.81*half/gmax
290 DO i = 1,nel
291 amax = amax1*max(cii(i,1),cii(i,2),cii(i,3))
292 eti(i) = max(one,amax)
293 gtmax(i) = gmax*eti(i) + gvmax
294 et(i) = eti(i)+etv
295 ENDDO
296c
297 !=======================================================================
298 ! - NEW STRESS TENSOR COMPUTATION
299 !=======================================================================
300c
301 ! Isochoric strain energy density derivation
302 ! Note: here, the table DWDL(:,J) corresponds to the product EVM(J)*DWDEVM(:,J)
303 DO j=1,3
304 dwdl(1:nel,j)=zero
305 DO ii=1,norder
306 DO i=1,nel
307 IF (evm(i,j)/=zero) THEN
308 IF (al(ii)/=zero) THEN
309 dwdl(i,j) = dwdl(i,j) + mu(ii) * exp(al(ii)*log(evm(i,j)))
310 ELSE
311 dwdl(i,j) = dwdl(i,j) + mu(ii)
312 ENDIF
313 ENDIF
314 ENDDO
315 ENDDO
316 ENDDO
317c
318 ! Volumic strain energy density derivation + trace of isochoric derivative
319 ! (Volumic strain energy density equation depending on formulation flag)
320 DO i=1,nel
321 ! -> Standard strain energy density
322 IF (iform == 1) THEN
323 dwdrv(i) = p(i)*(rv(i)- one)
324 ! -> Modified strain energy density
325 ELSEIF (iform == 2) THEN
326 dwdrv(i) = p(i)*(one - one/rv(i))
327 ENDIF
328 sumdwdl(i)=(dwdl(i,1)+dwdl(i,2)+dwdl(i,3))*third
329 ENDDO
330c
331 ! Cauchy principal stresses
332 DO j=1,3
333 DO i=1,nel
334 inve = one/rv(i)
335 sigprv(i,j) = (dwdl(i,j)-(sumdwdl(i)-rv(i)*dwdrv(i)))*inve
336 ENDDO
337 ENDDO
338c
339 ! Biot stress tensor for tension cutoff only
340 DO i=1,nel
341 t1(i,1) = rv(i)*sigprv(i,1)/ev(i,1)
342 t1(i,2) = rv(i)*sigprv(i,2)/ev(i,2)
343 t1(i,3) = rv(i)*sigprv(i,3)/ev(i,3)
344 ENDDO
345c
346 ! Tension cutoff stress
347 DO j=1,3
348 DO i=1,nel
349 IF(off(i)==zero.OR.t1(i,j)>abs(tensioncut))THEN
350 sigprv(i,1)=zero
351 sigprv(i,2)=zero
352 sigprv(i,3)=zero
353 off(i)=zero
354 ENDIF
355 ENDDO
356 ENDDO
357c
358 ! New stress tensor, soundspeed and user-variables
359 DO i=1,nel
360c
361 ! Transform principal Cauchy stresses to Global directions
362 signxx(i) = dirprv(i,1,1)*dirprv(i,1,1)*sigprv(i,1)
363 . + dirprv(i,1,2)*dirprv(i,1,2)*sigprv(i,2)
364 . + dirprv(i,1,3)*dirprv(i,1,3)*sigprv(i,3)
365c
366 signyy(i) = dirprv(i,2,2)*dirprv(i,2,2)*sigprv(i,2)
367 . + dirprv(i,2,3)*dirprv(i,2,3)*sigprv(i,3)
368 . + dirprv(i,2,1)*dirprv(i,2,1)*sigprv(i,1)
369c
370 signzz(i) = dirprv(i,3,3)*dirprv(i,3,3)*sigprv(i,3)
371 . + dirprv(i,3,1)*dirprv(i,3,1)*sigprv(i,1)
372 . + dirprv(i,3,2)*dirprv(i,3,2)*sigprv(i,2)
373c
374 signxy(i) = dirprv(i,1,1)*dirprv(i,2,1)*sigprv(i,1)
375 . + dirprv(i,1,2)*dirprv(i,2,2)*sigprv(i,2)
376 . + dirprv(i,1,3)*dirprv(i,2,3)*sigprv(i,3)
377c
378 signyz(i) = dirprv(i,2,2)*dirprv(i,3,2)*sigprv(i,2)
379 . + dirprv(i,2,3)*dirprv(i,3,3)*sigprv(i,3)
380 . + dirprv(i,2,1)*dirprv(i,3,1)*sigprv(i,1)
381c
382 signzx(i) = dirprv(i,3,3)*dirprv(i,1,3)*sigprv(i,3)
383 . + dirprv(i,3,1)*dirprv(i,1,1)*sigprv(i,1)
384 . + dirprv(i,3,2)*dirprv(i,1,2)*sigprv(i,2)
385c
386 ! Save user variables for post-processing
387 uvar(i,1) = max(sigprv(i,1),sigprv(i,2),sigprv(i,3))
388 uvar(i,2) = min(sigprv(i,1),sigprv(i,2),sigprv(i,3))
389 uvar(i,3) = off(i)
390 ! Soundspeed computation
391 ! Note: for Iform = 2, the constant bulk is assumed for the soundspeed computation
392 ! (if stability issue is observed, switch to non-linear bulk)
393 soundsp(i) = sqrt((two_third*gtmax(i)+p(i))/rho(i))
394 ! Viscosity parameter set to zero
395 viscmax(i) = zero
396 ENDDO
397c
398 !=======================================================================
399 ! - VISCOUS STRESSES COMPUTATION (PRONY SERIES)
400 !=======================================================================
401 IF (nprony > 0) THEN
402 IF (ismstr == 10 .or. ismstr == 12) THEN ! New formulation, using Dot(B)
403c---
404c F = deformation gradient
405c L = velocity gradient
406c L = D + W = sym(L) + skw(L) => LL = sym(L)
407c
408 DO i=1,nel
409 ll(i,1,1) = epsp1(i) ! EPSPXX
410 ll(i,2,2) = epsp2(i) ! EPSPYY
411 ll(i,3,3) = epsp3(i) ! EPSPZZ
412 ll(i,1,2) = epsp4(i) * half !(EPSPXY(I) + EPSPYX(I))/2
413 ll(i,2,3) = epsp5(i) * half !(EPSPYZ(I) + EPSPZY(I))/2
414 ll(i,1,3) = epsp6(i) * half !(EPSPZX(I) + EPSPXZ(I))/2
415 ll(i,2,1) = ll(i,1,2)
416 ll(i,3,1) = ll(i,1,3)
417 ll(i,3,2) = ll(i,2,3)
418c
419 f(i,1,1) = one + mfxx(i)
420 f(i,2,2) = one + mfyy(i)
421 f(i,3,3) = one + mfzz(i)
422 f(i,1,2) = mfxy(i)
423 f(i,2,3) = mfyz(i)
424 f(i,3,1) = mfzx(i)
425 f(i,2,1) = mfyx(i)
426 f(i,3,2) = mfzy(i)
427 f(i,1,3) = mfxz(i)
428 ENDDO
429c
430 CALL prodaat(f,bb,nel) ! B = F * FT
431c
432c Dev(B) = B * J^(-2/3), J = det(F)
433c
434 DO i=1,nel
435 bb(i,1,1) = bb(i,1,1) * j2third(i)
436 bb(i,2,2) = bb(i,2,2) * j2third(i)
437 bb(i,3,3) = bb(i,3,3) * j2third(i)
438 bb(i,1,2) = bb(i,1,2) * j2third(i)
439 bb(i,2,3) = bb(i,2,3) * j2third(i)
440 bb(i,1,3) = bb(i,1,3) * j2third(i)
441 bb(i,2,1) = bb(i,1,2)
442 bb(i,3,2) = bb(i,2,3)
443 bb(i,3,1) = bb(i,1,3)
444 ENDDO
445c
446 CALL prodmat(ll,bb,lb,nel) ! LB = L * B
447c
448 CALL prodmat(bb,ll,blt,nel) ! BLT = B * LT
449c
450 DO i=1,nel
451 dotb(i,1) = lb(i,1,1) + blt(i,1,1) ! xx
452 dotb(i,2) = lb(i,2,2) + blt(i,2,2) ! yy
453 dotb(i,3) = lb(i,3,3) + blt(i,3,3) ! zz
454 dotb(i,4) = lb(i,1,2) + blt(i,1,2) ! xy = yx
455 dotb(i,5) = lb(i,2,3) + blt(i,2,3) ! yz = zy
456 dotb(i,6) = lb(i,1,3) + blt(i,1,3) ! xz = zx
457 ENDDO
458c
459 DO j=1,6
460 svisc(1:nel,j) = zero
461 DO ii= 1,nprony
462 fac= -timestep/taux(ii)
463 DO i=1,nel
464 h0(i,j,ii) = uvar(i, 6 + (ii - 1)*6 + j)
465 h(i,j,ii) = exp(fac)*h0(i,j,ii) +
466 . exp(half*fac)*dotb(i,j)*timestep
467 uvar(i, 6 + (ii - 1)*6 + j)= h(i,j,ii)
468 ENDDO
469 ENDDO
470c
471c Kirchoff viscous stress
472c
473 DO ii = 1,nprony
474 DO i=1,nel
475 svisc(i,j) = svisc(i,j) + gi(ii)*h(i,j,ii)
476 ENDDO
477 ENDDO
478 ENDDO
479c
480c Kirchoff -> Cauchy visc stress : sig = T / J
481c
482 DO i=1,nel
483 inve = one/rv(i)
484 svisc(i,1) = svisc(i,1)*inve
485 svisc(i,2) = svisc(i,2)*inve
486 svisc(i,3) = svisc(i,3)*inve
487 svisc(i,4) = svisc(i,4)*inve
488 svisc(i,5) = svisc(i,5)*inve
489 svisc(i,6) = svisc(i,6)*inve
490c deviatoric part of visc stress
491 trace = (svisc(i,1) + svisc(i,2) + svisc(i,3)) * third
492 svisc(i,1) = svisc(i,1) - trace
493 svisc(i,2) = svisc(i,2) - trace
494 svisc(i,3) = svisc(i,3) - trace
495c total stress
496 signxx(i) = signxx(i) + svisc(i,1)
497 signyy(i) = signyy(i) + svisc(i,2)
498 signzz(i) = signzz(i) + svisc(i,3)
499 signxy(i) = signxy(i) + svisc(i,4)
500 signyz(i) = signyz(i) + svisc(i,5)
501 signzx(i) = signzx(i) + svisc(i,6)
502 ENDDO
503c
504 ELSE ! ISMSTR /= 10
505c
506 DO i=1,nel
507C
508 c0(i,1) = uvar(i,4)
509 c0(i,2) = uvar(i,5)
510 c0(i,3) = uvar(i,6)
511 c0(i,4) = uvar(i,7)
512 c0(i,5) = uvar(i,8)
513 c0(i,6) = uvar(i,9)
514C
515 cc = third*(evm(i,1)**2 + evm(i,2)**2 + evm(i,3)**2)
516 fft(1) = evm(i,1)**2 - cc
517 fft(2) = evm(i,2)**2 - cc
518 fft(3) = evm(i,3)**2 - cc
519c
520 c1(i,1) = dirprv(i,1,1)*dirprv(i,1,1)*fft(1)
521 . + dirprv(i,1,2)*dirprv(i,1,2)*fft(2)
522 . + dirprv(i,1,3)*dirprv(i,1,3)*fft(3)
523c
524 c1(i,2) = dirprv(i,2,2)*dirprv(i,2,2)*fft(2)
525 . + dirprv(i,2,3)*dirprv(i,2,3)*fft(3)
526 . + dirprv(i,2,1)*dirprv(i,2,1)*fft(1)
527c
528 c1(i,3) = dirprv(i,3,3)*dirprv(i,3,3)*fft(3)
529 . + dirprv(i,3,1)*dirprv(i,3,1)*fft(1)
530 . + dirprv(i,3,2)*dirprv(i,3,2)*fft(2)
531c
532 c1(i,4) = dirprv(i,1,1)*dirprv(i,2,1)*fft(1)
533 . + dirprv(i,1,2)*dirprv(i,2,2)*fft(2)
534 . + dirprv(i,1,3)*dirprv(i,2,3)*fft(3)
535c
536 c1(i,5) = dirprv(i,2,2)*dirprv(i,3,2)*fft(2)
537 . + dirprv(i,2,3)*dirprv(i,3,3)*fft(3)
538 . + dirprv(i,2,1)*dirprv(i,3,1)*fft(1)
539c
540 c1(i,6) = dirprv(i,3,3)*dirprv(i,1,3)*fft(3)
541 . + dirprv(i,3,1)*dirprv(i,1,1)*fft(1)
542 . + dirprv(i,3,2)*dirprv(i,1,2)*fft(2)
543C
544 uvar(i,4) = c1(i,1)
545 uvar(i,5) = c1(i,2)
546 uvar(i,6) = c1(i,3)
547 uvar(i,7) = c1(i,4)
548 uvar(i,8) = c1(i,5)
549 uvar(i,9) = c1(i,6)
550 ENDDO
551C
552 ! Viscous stresses computation
553 DO j=1,6
554 svisc(1:nel,j) = zero
555 DO ii= 1,nprony
556 fac= -timestep/taux(ii)
557 DO i=1,nel
558 h0(i,j,ii) = uvar(i, 12 + (ii - 1)*6 + j)
559 h(i,j,ii) = exp(fac)*h0(i,j,ii) +
560 . exp(half*fac)*(c1(i,j) - c0(i,j))
561 uvar(i, 12 + (ii - 1)*6 + j)= h(i,j,ii)
562 ENDDO
563 ENDDO
564C
565 DO ii = 1,nprony
566 DO i=1,nel
567 svisc(i,j) = svisc(i,j) + gi(ii)*h(i,j,ii)
568 ENDDO
569 ENDDO
570 ENDDO
571c
572 ! Add viscous stresses to the stress tensor
573 DO i=1,nel
574 inve = one/rv(i)
575 svisc(i,1) = svisc(i,1)*inve
576 svisc(i,2) = svisc(i,2)*inve
577 svisc(i,3) = svisc(i,3)*inve
578 svisc(i,4) = svisc(i,4)*inve
579 svisc(i,5) = svisc(i,5)*inve
580 svisc(i,6) = svisc(i,6)*inve
581C
582 signxx(i) = signxx(i) + svisc(i,1)
583 signyy(i) = signyy(i) + svisc(i,2)
584 signzz(i) = signzz(i) + svisc(i,3)
585 signxy(i) = signxy(i) + svisc(i,4)
586 signyz(i) = signyz(i) + svisc(i,5)
587 signzx(i) = signzx(i) + svisc(i,6)
588 ENDDO
589c
590 ENDIF ! ISMSTR
591 ENDIF
592C-----------
593 END
594
sigeps42
subroutine sigeps42(nel, nuparam, nuvar, nfunc, ifunc, npf, tf, time, timestep, uparam, rho0, rho, volume, eint, uvar, off, offg, soundsp, epsp1, epsp2, epsp3, epsp4, epsp5, epsp6, epsxx, epsyy, epszz, epsxy, epsyz, epszx, signxx, signyy, signzz, signxy, signyz, signzx, mfxx, mfxy, mfxz, mfyx, mfyy, mfyz, mfzx, mfzy, mfzz, viscmax, ismstr, et, ihet, epsth3, iexpan, niparam, iparam)
Definition sigeps42.F:45
min
#define min(a, b)
Definition macros.h:20
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
valpvecdp_v
subroutine valpvecdp_v(sig, val, vec, nel)
Definition sigeps33.F:2231
valpvec_v
subroutine valpvec_v(sig, val, vec, nel)
Definition sigeps33.F:1902