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

Go to the source code of this file.

Functions/Subroutines

subroutine sigeps42c (nel, nuparam, niparam, nuvar, ismstr, time, timestep, uparam, iparam, rho0, depsxx, depsyy, depsxy, depsyz, depszx, epsxx, epsyy, epsxy, thkn, thklyl, signxx, signyy, signxy, signyz, signzx, sigoyz, sigozx, soundsp, gs, uvar, off)

Function/Subroutine Documentation

◆ sigeps42c()

subroutine sigeps42c ( integer nel,
integer nuparam,
integer niparam,
integer nuvar,
integer ismstr,
time,
timestep,
uparam,
integer, dimension(niparam) iparam,
dimension(nel) rho0,
dimension(nel) depsxx,
dimension(nel) depsyy,
dimension(nel) depsxy,
dimension(nel) depsyz,
dimension(nel) depszx,
dimension(nel) epsxx,
dimension(nel) epsyy,
dimension(nel) epsxy,
dimension(nel) thkn,
dimension(nel) thklyl,
dimension(nel) signxx,
dimension(nel) signyy,
dimension(nel) signxy,
dimension(nel) signyz,
dimension(nel) signzx,
intent(in) sigoyz,
intent(in) sigozx,
dimension(nel) soundsp,
dimension(nel) gs,
uvar,
off )

Definition at line 28 of file sigeps42c.F.

36C-----------------------------------------------
37C I M P L I C I T T Y P E S
38C-----------------------------------------------
39#include "implicit_f.inc"
40C-----------------------------------------------
41C C O M M O N
42C-----------------------------------------------
43#include "param_c.inc"
44#include "com01_c.inc"
45C----------------------------------------------------------------
46C I N P U T A R G U M E N T S
47C----------------------------------------------------------------
48 INTEGER :: NEL
49 INTEGER :: NUPARAM
50 INTEGER :: NIPARAM
51 INTEGER :: NUVAR
52 INTEGER :: ISMSTR
53 INTEGER :: IPARAM(NIPARAM)
54 my_real :: uparam(nuparam)
55 my_real :: time,timestep
56 my_real ,DIMENSION(NEL) :: thkn,thklyl,rho0,gs,
57 . depsxx,depsyy,depsxy,depsyz,depszx,epsxx,epsyy,epsxy
58 my_real ,DIMENSION(NEL) ,INTENT(IN) :: sigoyz,sigozx
59C----------------------------------------------------------------
60C O U T P U T A R G U M E N T S
61C----------------------------------------------------------------
62 my_real ,DIMENSION(NEL) :: signxx,signyy,signxy,signyz,signzx,soundsp
63C----------------------------------------------------------------
64C I N P U T O U T P U T A R G U M E N T S
65C----------------------------------------------------------------
66 my_real :: uvar(nel,nuvar), off(nel)
67C----------------------------------------------------------------
68C L O C A L V A R I B L E S
69C----------------------------------------------------------------
70 INTEGER :: I,J,II,NPRONY,NORDER,ITER,JNV
71 my_real :: sum,fac,fscal,tenscut,sumdwdl,sumddwddl,partp
72 my_real :: e11,e22,e12,emax,gmax,gvmax,nu,rbulk,a11,pui11,pui22,alma1,alma2,alma3
73 my_real ,DIMENSION(3) :: ddwddl,dwdl
74 my_real ,DIMENSION(10) :: mu,al
75 my_real ,DIMENSION(100) :: gi,taux,h1,h2,h12,h10,h20,h120
76 my_real ,DIMENSION(NEL,100) :: h30,h31
77 my_real ,DIMENSION(NEL) :: rvt,c30,c31,dc3ev3,cd10,cd20,cd120
78 my_real ,DIMENSION(NEL) :: cp1,cp2,cd1,cd2,cd12,ea,invrv,invv3,dezz,rv,trav,rootv
79 my_real ,DIMENSION(NEL) :: s2,cs,kt3,epszz,rho
80 my_real ,DIMENSION(NEL,10) :: evma1,evma2,evma3,c11,c22,c12
81 my_real ,DIMENSION(NEL,3) :: t,sv,sigprv,evv,ev,evm
82 my_real ,DIMENSION(NEL,3,2) :: eigv(nel,3,2)
83C=======================================================================
84C SET INITIAL MATERIAL CONSTANTS
85
86 norder = iparam(1)
87 nprony = iparam(2)
88!
89 DO i=1,norder
90 mu(i) = uparam(i)
91 al(i) = uparam(i+10)
92 ENDDO
93 rbulk = uparam(21)
94 nu = uparam(22)
95 tenscut= uparam(23)
96 fscal = uparam(24)
97! ------------------
98 gmax = zero
99 gvmax = zero
100 DO i=1,nprony
101 gi(i) = uparam(24 + i)
102 taux(i) = uparam(24 + nprony + i)
103 gvmax = gvmax + gi(i)
104 ENDDO
105 DO i=1,norder
106 gmax = gmax + mu(i)*al(i)
107 ENDDO
108 gmax = gmax + gvmax
109C
110C User variables initialisation
111 IF (time == zero .AND. isigi == 0) THEN
112 DO i=1,nel
113 DO j=1,nuvar
114 uvar(i,j) = zero
115 ENDDO
116 uvar(i,3) = one
117 ENDDO
118 ELSEIF (time == zero ) THEN
119 DO i=1,nel
120 uvar(i,3) = one
121 ENDDO
122 ENDIF
123C principal stretch (def gradient eigenvalues)
124 DO i=1,nel
125 trav(i) = epsxx(i)+epsyy(i)
126 rootv(i) = sqrt((epsxx(i)-epsyy(i))*(epsxx(i)-epsyy(i))
127 . + epsxy(i)*epsxy(i))
128 evv(i,1) = half*(trav(i)+rootv(i))
129 evv(i,2) = half*(trav(i)-rootv(i))
130 evv(i,3) = zero
131 ea(i) = zero
132 ENDDO
133C-- avoid NaN---------
134 IF (ismstr == 10) THEN
135 DO i=1,nel
136 IF (min(evv(i,1),evv(i,2)) <= -one) THEN
137 evv(i,1) = zero
138 evv(i,2) = zero
139 off(i) = four_over_5
140 END IF
141 ENDDO
142 END IF
143C rot matrix (eigenvectors)
144 DO i=1,nel
145 IF (abs(evv(i,2)-evv(i,1)) < em10) THEN
146 eigv(i,1,1)=one
147 eigv(i,2,1)=one
148 eigv(i,3,1)=zero
149 eigv(i,1,2)=zero
150 eigv(i,2,2)=zero
151 eigv(i,3,2)=zero
152 ELSE
153 invrv(i) = one / rootv(i)
154 eigv(i,1,1) = (epsxx(i)-evv(i,2)) * invrv(i)
155 eigv(i,2,1) = (epsyy(i)-evv(i,2)) * invrv(i)
156 eigv(i,3,1) = (half*epsxy(i)) * invrv(i)
157 eigv(i,1,2) = (evv(i,1)-epsxx(i)) * invrv(i)
158 eigv(i,2,2) = (evv(i,1)-epsyy(i)) * invrv(i)
159 eigv(i,3,2) =-(half*epsxy(i)) * invrv(i)
160 ENDIF
161 ENDDO
162C Strain definition
163 IF (ismstr == 1 .OR. ismstr == 3 .OR. ismstr == 11) THEN ! engineering strain
164 DO i=1,nel
165 ev(i,1)=evv(i,1)+ one
166 ev(i,2)=evv(i,2)+ one
167 ev(i,3)=one/ev(i,1)/ev(i,2)
168 ENDDO
169 ELSEIF(ismstr == 10) THEN
170 DO i=1,nel
171 ev(i,1)=sqrt(evv(i,1)+ one)
172 ev(i,2)=sqrt(evv(i,2)+ one)
173 ev(i,3)=one/ev(i,1)/ev(i,2)
174 ENDDO
175 ELSE ! true strain
176 DO i=1,nel
177 ev(i,1)=exp(evv(i,1))
178 ev(i,2)=exp(evv(i,2))
179 ev(i,3)=one/ev(i,1)/ev(i,2)
180 ENDDO
181 ENDIF
182 IF (nprony > 0) ev(1:nel,3) = uvar(1:nel,3)
183 DO i=1,nel
184 IF (off(i)==zero.OR.off(i)==four_over_5) ev(i,1:3)=one
185 ENDDO
186C--------------------------------------
187C Newton method => Find EV(3) : T3(EV(3)) = 0
188C--------------------------------------
189 IF (nprony > 0) THEN
190C--------------------------------------
191C Newton method => Find EV(3) : T3(EV(3)) = 0
192C--------------------------------------
193 DO iter = 1,5
194! -----------------------
195 DO i=1,nel
196 rv(i) = ev(i,1)*ev(i,2)*ev(i,3)
197c---- normalized stretch => unified compressible/uncompressible formulation
198 ! RVT = RV(I)**(-THIRD)
199 ! RV = 0 --> RVT = 0
200 ! else --> RVT = exp( -third * ln( RV))
201 IF (rv(i)/= zero) THEN
202 rvt(i) = exp( (-third)* log(rv(i)) )
203 invrv(i) = one / rv(i)
204 ELSE
205 rvt(i) = zero
206 invrv(i) = zero
207 ENDIF
208 evm(i,1) = ev(i,1)*rvt(i)
209 evm(i,2) = ev(i,2)*rvt(i)
210 evm(i,3) = ev(i,3)*rvt(i)
211 ENDDO ! 1,NEL
212! -----------------------
213 ! partial derivatives of strain energy
214 DO j = 1,norder
215 DO i=1,nel
216 ! EVMA(I) = MU * EVM(I) ** AL(I) :
217 ! EVM = 0 --> EVMA(I) = 0
218 ! AL = 0 --> EVMA(I) = MU
219 ! else --> EVMA(I) = MU * exp(AL(I) * ln( EVM(I)))
220 IF(evm(i,1)==zero) THEN
221 evma1(i,j) = zero
222 ELSE
223 IF (al(j) == zero) THEN
224 evma1(i,j) = mu(j)
225 ELSE
226 evma1(i,j) = mu(j) * exp(al(j)* log(evm(i,1)) )
227 ENDIF
228 ENDIF
229!
230 IF(evm(i,2)==zero) THEN
231 evma2(i,j) = zero
232 ELSE
233 IF (al(j) == zero) THEN
234 evma2(i,j) = mu(j)
235 ELSE
236 evma2(i,j) = mu(j) * exp(al(j)* log(evm(i,2)) )
237 ENDIF
238 ENDIF
239!
240 IF(evm(i,3)==zero) THEN
241 evma3(i,j) = zero
242 ELSE
243 IF (al(j) == zero) THEN
244 evma3(i,j) = mu(j)
245 ELSE
246 evma3(i,j) = mu(j) * exp(al(j)* log(evm(i,3)) )
247 ENDIF
248 ENDIF
249 ENDDO ! 1,NEL
250 ENDDO ! J=1,NORDER
251! -----------------------
252 DO i=1,nel
253 dwdl(1) = zero
254 dwdl(2) = zero
255 dwdl(3) = zero
256 DO j=1,norder
257 dwdl(1) = dwdl(1) + evma1(i,j)
258 dwdl(2) = dwdl(2) + evma2(i,j)
259 dwdl(3) = dwdl(3) + evma3(i,j)
260 END DO
261 sumdwdl = (dwdl(1) + dwdl(2) + dwdl(3)) * third
262 partp = rbulk*(rv(i)- one)
263c---------
264c principal cauchy stress
265 IF (ev(i,3) == zero) THEN
266 invv3(i) = zero
267 ELSE
268 invv3(i) = one / ev(i,3)
269 ENDIF
270 t(i,1) = (dwdl(1) - sumdwdl) *invrv(i) + partp
271 t(i,2) = (dwdl(2) - sumdwdl) *invrv(i) + partp
272 t(i,3) = (dwdl(3) - sumdwdl) *invrv(i) + partp
273c---------
274 kt3(i) = -third*(dwdl(1) + dwdl(2)) + two_third*(dwdl(3))
275 kt3(i) = -ev(i,1)*ev(i,2)*kt3(i)*invrv(i)*invrv(i) + rbulk*ev(i,1)*ev(i,2)
276!
277 alma1 = zero
278 alma2 = zero
279 alma3 = zero
280 DO j=1,norder
281 alma1 = alma1 + al(j)*evma1(i,j)
282 alma2 = alma2 + al(j)*evma2(i,j)
283 alma3 = alma3 + al(j)*evma3(i,j)
284 END DO
285 kt3(i) = kt3(i) + (one_over_9*(alma1 + alma2 + four*(alma3))) * invrv(i)*invv3(i)
286C viscosty traitement
287 c30(i) = uvar(i,5)
288C
289 sum = third*(evm(i,1)**2 + evm(i,2)**2 + evm(i,3)**2)
290 c31(i) = evm(i,3)**2 - sum
291 ! SV3 = ZERO
292 dc3ev3(i) = four_over_3*rvt(i)*evm(i,3)-two_third*invv3(i)*
293 . (two_third*evm(i,3)**2 - third*evm(i,1)**2 - third*evm(i,2)**2)
294 ENDDO ! 1,NEL
295! -----------------------
296 jnv = 8
297 DO ii= 1,nprony
298 fac= -timestep/taux(ii)
299 h30(1:nel,ii) = uvar(1:nel,jnv + ii)
300 h31(1:nel,ii) = exp(fac)*h30(1:nel,ii)+ exp(half*fac)*(c31(1:nel) - c30(1:nel))
301 ENDDO
302C Kirchoff visco stress --->
303C PK2 stress, PK2 = F**(-1)*Taux* F**(-T)n
304C cauchy =Taux/RV is used here
305 DO ii = 1,nprony
306 fac= -timestep/taux(ii)
307 t(1:nel,3) = t(1:nel,3) + gi(ii)*h31(1:nel,ii)*invrv(1:nel)
308 kt3(1:nel) = kt3(1:nel) - gi(ii)*h31(1:nel,ii)*invv3(1:nel)*invrv(1:nel)
309 . + dc3ev3(1:nel)*gi(ii)*exp(half*fac)*invrv(1:nel)
310 ENDDO
311 DO i=1,nel
312 IF (off(i)==zero.OR.off(i)==four_over_5) cycle
313 IF (abs(kt3(i))>em20) ev(i,3) = ev(i,3) - t(i,3)/kt3(i)
314 ENDDO
315 ENDDO ! iteration
316! -----------------------
317C stored converged solution
318 jnv = 8
319 DO ii= 1,nprony
320 uvar(1:nel,jnv + ii) = h31(1:nel,ii)
321 ENDDO
322 DO i=1,nel
323 uvar(i,5) = c31(i)
324C
325 rv(i) = ev(i,1)*ev(i,2)*ev(i,3)
326! compute viscos stress
327 ! RVT = RV(I) ** (-THIRD) :
328 ! RV(I) = 0 --> RVT = 0
329 ! else --> RVT = exp((-THIRD) * ln(RV(I))
330 IF (rv(i) /= zero) THEN
331 rvt(i) = exp( (-third)* log(rv(i)) )
332 ELSE
333 rvt(i) = zero
334 ENDIF
335 evm(i,1) = ev(i,1)*rvt(i)
336 evm(i,2) = ev(i,2)*rvt(i)
337 evm(i,3) = ev(i,3)*rvt(i)
338C
339 cd10(i) = uvar(i,6)
340 cd20(i) = uvar(i,7)
341 cd120(i) = uvar(i,8)
342C
343 sum = third*(evm(i,1)**2 + evm(i,2)**2 + evm(i,3)**2)
344 cp1(i) = evm(i,1)**2 - sum
345 cp2(i) = evm(i,2)**2 - sum
346 cd1(i) = eigv(i,1,1)*cp1(i) + eigv(i,1,2)*cp2(i)
347 cd2(i) = eigv(i,2,1)*cp1(i) + eigv(i,2,2)*cp2(i)
348 cd12(i) = eigv(i,3,1)*cp1(i) + eigv(i,3,2)*cp2(i)
349 uvar(i,6) = cd1(i)
350 uvar(i,7) = cd2(i)
351 uvar(i,8) = cd12(i)
352 sv(i,1) = zero
353 sv(i,2) = zero
354 sv(i,3) = zero
355 ENDDO
356 jnv = 8 + nprony
357 DO ii= 1,nprony
358 DO i=1,nel
359 fac= -timestep/taux(ii)
360 h10(ii) = uvar(i,jnv + ii )
361 h20(ii) = uvar(i,jnv + nprony + ii )
362 h120(ii) = uvar(i,jnv + 2*nprony + ii )
363 h1(ii) = exp(fac)*h10(ii)+ exp(half*fac)*(cd1(i) - cd10(i))
364 h2(ii) = exp(fac)*h20(ii)+ exp(half*fac)*(cd2(i) - cd20(i))
365 h12(ii) = exp(fac)*h120(ii)+ exp(half*fac)*(cd12(i) - cd120(i))
366 uvar(i,jnv + ii )= h1(ii)
367 uvar(i,jnv + nprony + ii )= h2(ii)
368 uvar(i,jnv + 2*nprony + ii )= h12(ii)
369C Kirchoff visco stress
370 sv(i,1) = sv(i,1) + gi(ii)*h1(ii)
371 sv(i,2) = sv(i,2) + gi(ii)*h2(ii)
372 sv(i,3) = sv(i,3) + gi(ii)*h12(ii)
373 ENDDO ! 1,NEL
374 ENDDO ! NPRONY
375
376! -----------------------
377 ELSE ! lam3=1/lam1/lam2 (incompressible formulation) with out viscosity
378! -----------------------
379 DO j = 1,norder
380c---- normalized stretch => unified compressible/uncompressible formulation
381C---- partial derivatives of strain energy
382 ! EVMA(I) = MU * EV(I) ** AL(I) :
383 ! EV = 0 --> EVMA(I) = 0
384 ! AL = 0 --> EVMA(I) = MU
385 ! else --> EVMA(I) = MU * PUI
386 ! with PUI = exp(AL(I) * ln( EV(I)))
387 ! EVMA3 = MU * ( EV(1) * EV(2) )** (-AL(I))
388 ! --> EVMA3 = MU/ (PUI11 * PUI22)
389 DO i=1,nel
390 IF(ev(i,1)==zero) THEN
391 evma1(i,j) = zero
392 ELSE
393 IF (al(j) == zero) THEN
394 evma1(i,j) = mu(j)
395 pui11 = one
396 ELSE
397 pui11 = exp(al(j)* log(ev(i,1)) )
398 evma1(i,j) = mu(j) * pui11
399 ENDIF
400 ENDIF
401!
402 IF(ev(i,2)==zero) THEN
403 evma2(i,j) = zero
404 ELSE
405 IF (al(j) == zero) THEN
406 evma2(i,j) = mu(j)
407 pui22 = one
408 ELSE
409 pui22 = exp(al(j)* log(ev(i,2)) )
410 evma2(i,j) = mu(j) * pui22
411 ENDIF
412 ENDIF
413!
414 IF((ev(i,1)*ev(i,2))==zero) THEN
415 evma3(i,j) = zero
416 ELSE
417 IF (al(j) == zero) THEN
418 evma3(i,j) = mu(j)
419 ELSE
420 evma3(i,j) = mu(j)/(pui11*pui22)
421 ENDIF
422 ENDIF
423 ENDDO ! 1,NEL
424 ENDDO ! 1,NORDER
425! -----------------------
426 DO i=1,nel
427 dwdl(1) = zero
428 dwdl(2) = zero
429 dwdl(3) = zero
430 DO j=1,norder
431 dwdl(1) = dwdl(1) + evma1(i,j)
432 dwdl(2) = dwdl(2) + evma2(i,j)
433 dwdl(3) = dwdl(3) + evma3(i,j)
434 END DO
435c---------
436c principal cauchy stress
437 t(i,1) = dwdl(1) - dwdl(3)
438 t(i,2) = dwdl(2) - dwdl(3)
439 t(i,3) = zero
440 ENDDO
441! -----------------------
442! compute rigidity
443 DO j = 1,norder
444 ! CXX = 1/2 * MU * [1/2 * AL - 1] * EVX ** (AL - 4)
445 ! CXY = 1/2 * MU * [1/2 * AL + 1] * (EVX * EVY) **(-AL - 4)
446 ! AL - 4 = 0 --> CXX = 1/2 * MU * [1/2 * AL - 1]
447 ! --> CXY = 1/2 * MU * [1/2 * AL + 1] * (EVX * EVY) **(-8)
448 ! AL = 0 --> CXX = 1/2 * MU * [1/2 * AL - 1]* EVX ** (- 4)
449 ! --> CXY = 1/2 * MU * [1/2 * AL + 1] * (EVX * EVY) **(-4)
450 ! else --> CXX = 1/2 * MU * [1/2 * AL - 1]* PUIXX / EVX**(4)
451 ! with PUIXX = exp( AL * ln(EVX) )
452 ! --> CXY = 1/2 * MU * [1/2 * AL + 1] / ( PUIXX*PUIYY*(EVX**(4)*EVY**(4) )
453 DO i=1,nel
454 IF (al(j) == four) THEN
455 c11(i,j) = half*mu(j)*(half*al(j)-one)
456 c22(i,j) = half*mu(j)*(half*al(j)-one)
457 c12(i,j) = half*mu(j)*(half*al(j) + one) / (ev(i,1)*ev(i,2))**(8)
458 ELSEIF (al(j) == zero) THEN
459 c11(i,j) = half*mu(j)*(half*al(j)-one) / ev(i,1)**(4)
460 c22(i,j) = half*mu(j)*(half*al(j)-one) / ev(i,2)**(4)
461 c12(i,j) = half*mu(j)*(half*al(j) + one) / (ev(i,1)*ev(i,2))**(4)
462 ELSE
463 IF(ev(i,1) /= zero) THEN
464 pui11 = exp((al(j) )* log(ev(i,1)) )
465 c11(i,j) = half*mu(j)*(half*al(j) - one)*pui11 / ev(i,1)**(4)
466 ELSE
467 c11(i,j) = zero
468 pui11 = one
469 ENDIF
470 IF(ev(i,2) /= zero) THEN
471 pui22 = exp((al(j))* log(ev(i,2)) )
472 c22(i,j) = half*mu(j)*(half*al(j) - one)*pui22 / ev(i,2)**(4)
473 ELSE
474 c22(i,j) = zero
475 pui22 = one
476 ENDIF
477 IF(ev(i,1)*ev(i,2) /= zero) THEN
478 c12(i,j) = half*mu(j)*(half*al(j) + one) /
479 . ( (pui11 * pui22)*(ev(i,1)**4 * ev(i,2)**4) )
480 ELSE
481 c12(i,j) = zero
482 ENDIF
483 ENDIF
484 ENDDO ! 1,NEL
485 ENDDO ! 1,NORDER
486! -----------------------
487 DO i=1,nel
488 e11 = zero
489 e22 = zero
490 e12 = zero
491 DO j=1,norder
492 e11 = e11 + c11(i,j)
493 e22 = e22 + c22(i,j)
494 e12 = e12 + c12(i,j)
495 END DO
496 e11 = e11 + e12 * ev(i,2)**4
497 e22 = e22 + e12 * ev(i,1)**4
498 ea(i) = max(e11,e22)
499 sv(i,1) = zero
500 sv(i,2) = zero
501 sv(i,3) = zero
502 ENDDO ! 1,NEL
503! -----------------------
504 ENDIF ! PRONY
505C-------------------------------------------------------------
506c tension cut
507 DO i=1,nel
508 IF (off(i) /= zero .AND.
509 . (t(i,1) > abs(tenscut) .OR. t(i,2) > abs(tenscut))) THEN
510 t(i,1) = zero
511 t(i,2) = zero
512 t(i,3) = zero
513 off(i) = four_over_5
514 ENDIF
515 ENDDO
516C-------------------------------------------------------------
517C transform principal Cauchy stress to global directions
518C-------------------------------------------------------------
519 IF (ismstr == 1 .OR. ismstr == 3 .OR. ismstr == 11) THEN ! engineering strain
520 DO i=1,nel
521 epszz(i) = ev(i,3) - one
522 uvar(i,3) = ev(i,3)
523 ENDDO
524 ELSEIF (ismstr == 10) THEN ! left gauchy-green strain
525 DO i=1,nel
526 epszz(i) = ev(i,3) - one
527 uvar(i,3) = ev(i,3)
528 ENDDO
529 ELSE ! true strain
530 DO i=1,nel
531 epszz(i) = log(ev(i,3))
532 uvar(i,3) = ev(i,3)
533 ENDDO
534 ENDIF
535c
536 DO i=1,nel
537 rv(i) = ev(i,1)*ev(i,2)*ev(i,3)
538 IF (rv(i) /= zero) THEN
539 invrv(i) = one / rv(i)
540 ELSE
541 invrv(i) = zero
542 ENDIF
543c
544 dezz(i) = -nu/(one-nu)*(depsxx(i)+depsyy(i))
545!! DEZZ(I) = EPSZZ(I) - UVAR(I,4)
546 signxx(i) = eigv(i,1,1)*t(i,1) + eigv(i,1,2)*t(i,2) + sv(i,1)*invrv(i)
547 signyy(i) = eigv(i,2,1)*t(i,1) + eigv(i,2,2)*t(i,2) + sv(i,2)*invrv(i)
548 signxy(i) = eigv(i,3,1)*t(i,1) + eigv(i,3,2)*t(i,2) + sv(i,3)*invrv(i)
549 signyz(i) = sigoyz(i)+gs(i)*depsyz(i)
550 signzx(i) = sigozx(i)+gs(i)*depszx(i)
551 rho(i) = rho0(i)*invrv(i)
552 thkn(i) = thkn(i) + dezz(i)*thklyl(i)*off(i)
553 uvar(i,4) = epszz(i)
554C
555 emax = gmax*(one + nu)
556 emax = max(emax,ea(i))
557 a11 = emax/(one - nu**2)
558 soundsp(i)= sqrt(a11/rho0(i))
559 ENDDO
560C-----------
561 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21