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

Go to the source code of this file.

Functions/Subroutines

subroutine sigeps62c (nel, nuparam, nuvar, nfunc, ifunc, npf, npt0, ilayer, iflag, tf, time, timestep, uparam, rho0, area, eint, thklyl, epspxx, epspyy, epspxy, epspyz, epspzx, depsxx, depsyy, depsxy, depsyz, depszx, epsxx, epsyy, epsxy, epsyz, epszx, sigoxx, sigoyy, sigoxy, sigoyz, sigozx, signxx, signyy, signxy, signyz, signzx, sigvxx, sigvyy, sigvxy, sigvyz, sigvzx, soundsp, viscmax, thkn, uvar, off, ngl, ismstr, ipm, gs)

Function/Subroutine Documentation

◆ sigeps62c()

subroutine sigeps62c ( integer nel,
integer nuparam,
integer nuvar,
integer nfunc,
integer, dimension(nfunc) ifunc,
integer, dimension(*) npf,
integer npt0,
integer ilayer,
integer, dimension(*) iflag,
tf,
time,
timestep,
uparam,
rho0,
area,
eint,
thklyl,
epspxx,
epspyy,
epspxy,
epspyz,
epspzx,
depsxx,
depsyy,
depsxy,
depsyz,
depszx,
epsxx,
epsyy,
epsxy,
epsyz,
epszx,
sigoxx,
sigoyy,
sigoxy,
sigoyz,
sigozx,
signxx,
signyy,
signxy,
signyz,
signzx,
sigvxx,
sigvyy,
sigvxy,
sigvyz,
sigvzx,
soundsp,
viscmax,
thkn,
uvar,
off,
integer, dimension(nel) ngl,
integer ismstr,
integer, dimension(npropmi,*) ipm,
gs )

Definition at line 30 of file sigeps62c.F.

43C-----------------------------------------------
44C I M P L I C I T T Y P E S
45C-----------------------------------------------
46#include "implicit_f.inc"
47C-----------------------------------------------
48C C O M M O N
49C-----------------------------------------------
50#include "param_c.inc"
51C----------------------------------------------------------------
52C I N P U T A R G U M E N T S
53C----------------------------------------------------------------
54 INTEGER NEL,NUPARAM, NUVAR,ISMSTR, IPM(NPROPMI,*),MAT(NEL),
55 . NPT0,ILAYER, IFLAG(*),NGL(NEL)
56 my_real
57 . time,timestep,uparam(nuparam),thkn(nel),thklyl(nel),
58 . rho0(nel),area(nel),eint(nel,2),gs(nel),
59 . epspxx(nel),epspyy(nel),epspxy(nel),epspyz(nel),epspzx(nel),
60 . depsxx(nel),depsyy(nel),depsxy(nel),depsyz(nel),depszx(nel),
61 . epsxx(nel),epsyy(nel),epsxy(nel),epsyz(nel),epszx(nel),
62 . sigoxx(nel),sigoyy(nel),sigoxy(nel),sigoyz(nel),sigozx(nel)
63C----------------------------------------------------------------
64C O U T P U T A R G U M E N T S
65C----------------------------------------------------------------
66 my_real
67 . signxx(nel),signyy(nel),signxy(nel),signyz(nel),signzx(nel),
68 . sigvxx(nel),sigvyy(nel),sigvxy(nel),sigvyz(nel),sigvzx(nel),
69 . soundsp(nel),viscmax(nel),et(nel)
70C----------------------------------------------------------------
71C I N P U T O U T P U T A R G U M E N T S
72C----------------------------------------------------------------
73 my_real
74 . uvar(nel,nuvar), off(nel)
75C----------------------------------------------------------------
76C VARIABLES FOR FUNCTION INTERPOLATION
77C----------------------------------------------------------------
78 INTEGER NPF(*), NFUNC, IFUNC(NFUNC)
79 my_real finter,fintte,tf(*),fint2v
80 EXTERNAL finter,fintte
81C----------------------------------------------------------------
82C L O C A L V A R I B L E S
83C----------------------------------------------------------------
84 INTEGER I,J,K,NORDRE,NPRONY,ITER,IVISC,JNV
85 my_real
86 . nu,gmax,rbulk,beta(100),gamainf,dvisc,sum,fac,rvt,evl,
87 . ev1,ev2,ev3,kt3(nel),invr(nel),dsdev(nel),dsvol(nel),
88 . dq1(nel),dq2(nel),dq3(nel),sdh01(nel),
89 . sdh02(nel),sdh03(nel),sdh1(nel),sdh2(nel),sdh3(nel),
90 . dqdev(nel),dqidev,inv31,inv32(nel),inv33(nel),
91 . rvd(nel),h10,h20,h30,h1,h2,h3,hp1,hp2,hp3,hd1,hd2,hd3,a11,
92 . emax,sdhp1,sdhp2,inv11,inv21
93 my_real
94 . evv(nel,3),ev(nel,3),evm(nel),dwdl(3),ec(nel,3),
95 . rho(nel),gama(nel),pres(nel),rv(nel),p,sp(nel,3),
96 . t(nel,3),sv(nel,3),s(nel,3),s0(nel,3),sd(nel,3),sd0(nel,3),
97 . trav(nel),rootv(nel),eigv(nel,3,2),dezz(nel),
98 . ev1_sav(nel,100),ev2_sav(nel,100)
99 my_real
100 . mu(100),al(100),d(100),
101 . taux(100)
102C=======================================================================
103C SET INITIAL MATERIAL CONSTANTS
104 nu = uparam(1)
105 nordre = nint(uparam(2))
106 nprony = nint(uparam(3))
107 gamainf = uparam(4)
108 rbulk = uparam(5)
109 dvisc = uparam(6)
110 gmax = zero
111 DO i= 1,nordre
112 mu(i) = uparam(6 + i)
113 al(i) = uparam(6 + nordre + i)
114 beta(i) = uparam(2*nordre + 2*nprony + 7 + i)
115 gmax = gmax + mu(i)
116 ENDDO
117 ivisc = 0
118 IF (nprony /= zero) THEN
119 DO i=1,nprony
120 gama(i) = uparam(6 + 2*nordre + i)
121 taux(i) = uparam(6 + 2*nordre + nprony + i)
122 ENDDO
123 ivisc = 1
124 ENDIF
125 gmax = two*gmax
126C User variables initialisation
127 IF (time == zero) THEN
128 DO i = 1, nel
129 DO j = 1,nuvar
130 uvar(i,j) = zero
131 ENDDO
132 uvar(i,4) = one
133 ENDDO
134 ENDIF
135C principal stretch (def gradient eigenvalues)
136 DO i=1,nel
137 trav(i) = epsxx(i)+epsyy(i)
138 rootv(i) = sqrt((epsxx(i)-epsyy(i))*(epsxx(i)-epsyy(i))
139 . + epsxy(i)*epsxy(i))
140 evv(i,1) = half*(trav(i)+rootv(i))
141 evv(i,2) = half*(trav(i)-rootv(i))
142 evv(i,3) = zero
143 ENDDO
144C rot matrix (eigenvectors)
145 DO i=1,nel
146 IF (abs(evv(i,2)-evv(i,1)) < em10) THEN
147 eigv(i,1,1)=one
148 eigv(i,2,1)=one
149 eigv(i,3,1)=zero
150 eigv(i,1,2)=zero
151 eigv(i,2,2)=zero
152 eigv(i,3,2)=zero
153 ELSE
154 eigv(i,1,1) = (epsxx(i)-evv(i,2))/rootv(i)
155 eigv(i,2,1) = (epsyy(i)-evv(i,2))/rootv(i)
156 eigv(i,3,1) = (half*epsxy(i)) /rootv(i)
157 eigv(i,1,2) = (evv(i,1)-epsxx(i))/rootv(i)
158 eigv(i,2,2) = (evv(i,1)-epsyy(i))/rootv(i)
159 eigv(i,3,2) =-(half*epsxy(i)) /rootv(i)
160 ENDIF
161 ENDDO
162C Strain definition
163 IF (ismstr == 1 .OR. ismstr == 3) 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)=uvar(i,4)
168 ENDDO
169 ELSE ! true strain
170 DO i=1,nel
171 ev(i,1)=exp(evv(i,1))
172 ev(i,2)=exp(evv(i,2))
173 ev(i,3)=uvar(i,4)
174 ENDDO
175 ENDIF
176C--------------------------------------
177C Newton method => Find EV(3) : T3(EV(3)) = 0
178C--------------------------------------
179 IF(ivisc == 0) THEN
180 DO k=1,nordre
181 DO i=1,nel
182 ev1 = zero
183 IF(ev(i,1) > zero) THEN
184 ev1 = one
185 IF(al(k) /= zero)ev1 = exp(al(k)*log(ev(i,1)))
186 ENDIF
187 ev2 = zero
188 IF(ev(i,2)> zero) THEN
189 ev2 = one
190 IF(al(k) /= zero) ev2 = exp(al(k)*log(ev(i,2)))
191 ENDIF
192 ev1_sav(i,k) = ev1
193 ev2_sav(i,k) = ev2
194 ENDDO
195 ENDDO
196! -----------------------
197 DO iter = 1,5
198C
199! RVT = RV(I)**(-BETA*AL(K))
200! if rv > 0 --> pui = exp(-beta * ln(ev))
201! else pui = 0
202 rv(1:nel) = ev(1:nel,1)*ev(1:nel,2)*ev(1:nel,3)
203 t(1:nel,3) = zero
204 kt3(1:nel) = zero
205! -----------------------
206 DO k=1,nordre
207 DO i=1,nel
208 ev1 = ev1_sav(i,k)
209 ev2 = ev2_sav(i,k)
210 ev3 = zero
211 IF(ev(i,3)> zero) THEN
212 ev3 = one
213 IF(al(k) /= zero) ev3 = exp(al(k)*log(ev(i,3)))
214 ENDIF
215 rvt = zero
216 IF(rv(i) > zero) THEN
217 rvt = one
218 IF(al(k)/= zero) rvt = exp((-beta(k)*al(k))*log(rv(i)))
219 ENDIF
220C
221 fac = two*mu(k)/al(k)
222C---- cauchy stress
223 t(i,3) = t(i,3) + fac*(ev3 - rvt)/rv(i)
224
225 kt3(i) = kt3(i) +
226 . (-fac *(ev3 - rvt) +
227 . two*mu(k)*(ev3 + beta(k)*rvt))/ev(i,3)/rv(i)
228 ENDDO
229 ENDDO ! Ordre
230! -----------------------
231 DO i=1,nel
232 IF(kt3(i) /= zero) ev(i,3)= ev(i,3) - t(i,3) / kt3(i)
233 ENDDO
234 ENDDO ! Ite
235! -----------------------
236
237 uvar(1:nel,4) = ev(1:nel,3)
238C compute cauchy stress using the converged value
239 rv(1:nel) = ev(1:nel,1)*ev(1:nel,2)*ev(1:nel,3)
240C
241 t(1:nel,1) = zero
242 t(1:nel,2) = zero
243 t(1:nel,3) = zero
244! -----------------------
245 DO k=1,nordre
246 DO i=1,nel
247 ev1 = ev1_sav(i,k)
248 ev2 = ev2_sav(i,k)
249 ev3 = zero
250 IF(ev(i,3)> zero) THEN
251 ev3 = one
252 IF(al(k) /= zero) ev3 = exp(al(k)*log(ev(i,3)))
253 ENDIF
254 rvt = zero
255 IF(rv(i) > zero) THEN
256 rvt = one
257 IF(al(k)/= zero) rvt = exp((-beta(k)*al(k))*log(rv(i)))
258 ENDIF
259
260 fac = two*mu(k)/al(k)
261 t(i,1) = t(i,1) + fac*(ev1 - rvt)/rv(i)
262 t(i,2) = t(i,2) + fac*(ev2 - rvt)/rv(i)
263 t(i,3) = t(i,3) + fac*(ev3 - rvt)/rv(i)
264 ENDDO
265 ENDDO
266! -----------------------
267 ELSE ! with viscosity
268C-----------------------------------------------------------
269C Maxwell model -
270C-----------------------------------------------------------
271 DO k=1,nordre
272 DO i=1,nel
273 ev1 = zero
274 IF(ev(i,1) > zero) THEN
275 ev1 = one
276 IF(al(k) /= zero)ev1 = exp(al(k)*log(ev(i,1)))
277 ENDIF
278 ev2 = zero
279 IF(ev(i,2)> zero) THEN
280 ev2 = one
281 IF(al(k) /= zero) ev2 = exp(al(k)*log(ev(i,2)))
282 ENDIF
283 ev1_sav(i,k) = ev1
284 ev2_sav(i,k) = ev2
285 ENDDO
286 ENDDO
287! -----------------------
288 DO iter = 1,5
289 rv(1:nel) = ev(1:nel,1)*ev(1:nel,2)*ev(1:nel,3)
290 invr(1:nel) = one/max(em20,rv(1:nel))
291 dsdev(1:nel) = zero
292 dsvol(1:nel) = zero
293 kt3(1:nel) = zero
294 t(1:nel,1) = zero
295 t(1:nel,2) = zero
296 t(1:nel,3) = zero
297 dq1(1:nel) = zero
298 dq2(1:nel) = zero
299 dq3(1:nel) = zero
300 DO k=1,nordre
301 DO i=1,nel
302 ev1 = ev1_sav(i,k)
303 ev2 = ev2_sav(i,k)
304 ev3 = zero
305 IF(ev(i,3)> zero) THEN
306 ev3 = one
307 IF(al(k) /= zero) ev3 = exp(al(k)*log(ev(i,3)))
308 ENDIF
309C
310 rvt = zero
311 IF(rv(i) > zero) THEN
312 rvt = one
313 IF(al(k)/= zero) rvt = exp((-beta(k)*al(k))*log(rv(i)))
314 ENDIF
315C
316 fac = two*mu(k)/al(k)
317C---- cauchy stress
318 t(i,1) = t(i,1) + fac*(ev1 - rvt)*invr(i)
319 t(i,2) = t(i,2) + fac*(ev2 - rvt)*invr(i)
320 t(i,3) = t(i,3) + fac*(ev3 - rvt)*invr(i)
321C
322 dsvol(i) = dsvol(i) - two_third*fac*(ev1 + ev2 + ev3)
323 . + two*fac*rvt + two_third*mu(k)*ev3
324 . + two*beta(k)*mu(k)*rvt
325C
326 dsdev(i) = dsdev(i) - two*fac*(two_third*ev3 - third*(ev1 + ev2))
327 . + four_over_3*mu(k)*ev3
328 dq1(i) = dq1(i) - two_third*mu(k)*ev1
329 dq2(i) = dq2(i) - two_third*mu(k)*ev2
330 dq3(i) = dq3(i) + four_over_3*mu(k)*ev3
331 ENDDO
332 ENDDO ! Nordre
333
334
335 DO i=1,nel
336 inv31 = one/ev(i,3)
337 inv32(i) = inv31*inv31
338 inv33(i) = inv32(i)*inv31
339 inv11 = one/max(em20,ev(i,1))
340 inv21 = one/max(em20,ev(i,2))
341C
342 dsvol(i) = dsvol(i)*inv33(i)
343 dsdev(i) = dsdev(i)*inv33(i)
344 dq1(i) = dq1(i)*inv31*inv11*inv11
345 dq2(i) = dq2(i)*inv31*inv21*inv21
346 dq2(i) = dq3(i)*inv31
347C
348C Deviatoric stress + pressure
349C
350 p = third*(t(i,1) + t(i,2) + t(i,3))
351C PK2 pressure stress
352 ec(i,1) = ev(i,1)*ev(i,1)
353 ec(i,2) = ev(i,2)*ev(i,2)
354 ec(i,3) = ev(i,3)*ev(i,3)
355C PK2 pressure stress
356 sp(i,1)= rv(i)*p/max(em20,ec(i,1))
357 sp(i,2)= rv(i)*p/max(em20,ec(i,2))
358 sp(i,3)= rv(i)*p/max(em20,ec(i,3))
359
360C Deviatoric PK2 stress
361 sd(i,1)= rv(i)*(t(i,1) - p)/max(em20,ec(i,1))
362 sd(i,2)= rv(i)*(t(i,2) - p)/max(em20,ec(i,2))
363 sd(i,3)= rv(i)*(t(i,3) - p)/max(em20,ec(i,3))
364 rvd(i) = zero
365 IF(rv(i)>zero) rvd(i) = exp( two_third*log(rv(i)))
366C
367 sdhp1 = rvd(i)*sd(i,1)
368 sdhp2 = rvd(i)*sd(i,2)
369 sdh3(i) = rvd(i)*sd(i,3)
370
371 dq1(i) = dq1(i)*rvd(i) + two_third*sdhp1*inv31
372 dq2(i) = dq2(i)*rvd(i) + two_third*sdhp2*inv31
373 dq2(i) = dq2(i)*rvd(i) + two_third*sdh3(i)*ev(i,3)
374
375C old Deviatoric PK2 stress in the global frame
376 sdh01(i) = uvar(i,1)
377 sdh02(i) = uvar(i,2)
378 sdh03(i) = uvar(i,3)
379C
380 sdh1(i) = eigv(i,1,1)*sdhp1 + eigv(i,1,2)*sdhp2
381 sdh2(i) = eigv(i,2,1)*sdhp1 + eigv(i,2,2)*sdhp2
382 !SDH3(I) = SDH3(I)
383
384 dqdev(i) = zero
385 rvd(i) = zero
386 IF(rv(i) > zero)rvd(i) = exp( (-two_third)*log(rv(i)) )
387 s(i,3) = sp(i,3) + gamainf*sd(i,3)
388 ENDDO
389 jnv = 4
390C viscosity
391 DO k = 1,nprony
392 dqidev = zero
393 fac= -timestep/taux(k)
394 DO i=1,nel
395 h10 = uvar(i,jnv + k )
396 h20 = uvar(i,jnv + nprony + k )
397 h30 = uvar(i,jnv + 2*nprony + k )
398C
399 h1 = exp(fac)*h10+ exp(half*fac)*(sdh1(i) - sdh01(i))
400 h2 = exp(fac)*h20+ exp(half*fac)*(sdh2(i) - sdh02(i))
401 h3 = exp(fac)*h30+ exp(half*fac)*(sdh3(i) - sdh03(i))
402
403C H in the principal frame
404 hp1 = eigv(i,1,1)*h1 + eigv(i,2,1)*h2
405 hp2 = eigv(i,1,2)*h1 + eigv(i,2,2)*h2
406
407C H Deviatoric part
408 fac = third*(hp1*ec(i,1) + hp2*ec(i,2) + h3*ec(i,3))
409 hd1 = hp1 - fac/max(em20,ec(i,1))
410 hd2 = hp2 - fac/max(em20,ec(i,2))
411 hd3 = h3 - fac/max(em20,ec(i,3))
412C compute stress
413 s(i,3) = s(i,3) + gama(k)*rvd(i)*hd3
414 dqidev = exp(half*fac)*(dsdev(i)
415 . -third*(ec(i,1)*dq1(i)+ec(i,2)*dq2(i)+ec(i,3)*dq3(i))*inv32(i))
416 . +two_third*(ec(i,1)*hp1+ ec(i,2)*hp2 + h3*ec(i,3))*inv33(i)
417 dqdev(i)=dqdev(i) + rvd(i)*(-gama(k)*dqidev + two_third*gama(k)*hd3 )
418 ENDDO ! 1,NEL
419 ENDDO ! NPRONY
420 kt3(1:nel) = dsvol(1:nel) + gamainf*dsdev(1:nel) + dqdev(1:nel)
421C update of lam_3
422 DO i=1,nel
423 IF(kt3(i) /= zero) ev(i,3) = ev(i,3) - s(i,3)/ kt3(i)
424 ENDDO
425 ENDDO ! iter
426
427
428
429
430
431C stocked value
432 uvar(1:nel,4) = ev(1:nel,3)
433Compute stress after convergency
434 ! RVT = RV(I)**(-BETA*AL(K))
435! if rv > 0 --> pui = exp(-beta * ln(ev))
436! else pui = 0
437 rv(1:nel) = ev(1:nel,1)*ev(1:nel,2)*ev(1:nel,3)
438C
439 t(1:nel,1) = zero
440 t(1:nel,2) = zero
441 t(1:nel,3) = zero
442 DO k=1,nordre
443 DO i=1,nel
444 ev1 = ev1_sav(i,k)
445 ev2 = ev2_sav(i,k)
446C
447 rvt = zero
448 IF(rv(i) > zero) THEN
449 rvt = one
450 IF(al(k)/= zero) rvt = exp((-beta(k)*al(k))*log(rv(i)))
451 ENDIF
452C
453 ev3 = zero
454 IF(ev(i,3)> zero) THEN
455 ev3 = one
456 IF(al(k) /= zero) ev3 = exp(al(k)*log(ev(i,3)))
457 ENDIF
458
459 fac = two*mu(k)/al(k)
460C---- cauchy stress
461 t(i,1) = t(i,1) + fac*(ev1 - rvt)/rv(i)
462 t(i,2) = t(i,2) + fac*(ev2 - rvt)/rv(i)
463 t(i,3) = t(i,3) + fac*(ev3 - rvt)/rv(i)
464 ENDDO
465 ENDDO ! NORDRE
466
467 DO i=1,nel
468C
469C Deviatoric stress + pressure
470C
471 p = third*(t(i,1) + t(i,2) + t(i,3))
472C
473 ec(i,1) = ev(i,1)*ev(i,1)
474 ec(i,2) = ev(i,2)*ev(i,2)
475 ec(i,3) = ev(i,3)*ev(i,3)
476C PK2 pressure stress
477 sp(i,1)= rv(i)*p/max(em20,ec(i,1))
478 sp(i,2)= rv(i)*p/max(em20,ec(i,2))
479 sp(i,3)= rv(i)*p/max(em20,ec(i,3))
480C Deviatoric PK2 stress
481 sd(i,1)= rv(i)*(t(i,1) - p)/max(em20,ec(i,1))
482 sd(i,2)= rv(i)*(t(i,2) - p)/max(em20,ec(i,2))
483 sd(i,3)= rv(i)*(t(i,3) - p)/max(em20,ec(i,3))
484 rvd(i) = zero
485 IF(rv(i)>zero) rvd(i) = exp( two_third*log(rv(i)))
486C
487 sdhp1 = rvd(i)*sd(i,1)
488 sdhp2 = rvd(i)*sd(i,2)
489 sdh3(i) = rvd(i)*sd(i,3)
490C
491C compute cauchy stress using the converged value
492C
493
494C old Deviatoric PK2 stress in the global frame
495 sdh01(i) = uvar(i,1)
496 sdh02(i) = uvar(i,2)
497 sdh03(i) = uvar(i,3)
498C
499 sdh1(i) = eigv(i,1,1)*sdhp1 + eigv(i,1,2)*sdhp2
500 sdh2(i) = eigv(i,2,1)*sdhp1 + eigv(i,2,2)*sdhp2
501 !SDH3(I) = SDH3(I)
502
503 uvar(i,1) = sdh1(i)
504 uvar(i,2) = sdh2(i)
505 uvar(i,3) = sdh3(i)
506
507 rvd(i) = zero
508 IF(rv(i) > zero)rvd(i) = exp( (-two_third)*log(rv(i)) )
509 s(i,1) = sp(i,1) + gamainf*sd(i,1)
510 s(i,2) = sp(i,2) + gamainf*sd(i,2)
511 s(i,3) = sp(i,3) + gamainf*sd(i,3)
512 ENDDO
513C
514 jnv = 4
515 DO k = 1,nprony
516 DO i=1,nel
517 fac= -timestep/taux(k)
518 h10 = uvar(i,jnv + k )
519 h20 = uvar(i,jnv + nprony + k )
520 h30 = uvar(i,jnv + 2*nprony + k )
521 h1 = exp(fac)*h10+ exp(half*fac)*(sdh1(i) - sdh01(i))
522 h2 = exp(fac)*h20+ exp(half*fac)*(sdh2(i) - sdh02(i))
523 h3 = exp(fac)*h30+ exp(half*fac)*(sdh3(i) - sdh03(i))
524C
525 uvar(i,jnv + k )= h1
526 uvar(i,jnv + nprony + k )= h2
527 uvar(i,jnv + 2*nprony + k )= h3
528C H in the principal frame
529 hp1 = eigv(i,1,1)*h1 + eigv(i,2,1)*h2
530 hp2 = eigv(i,1,2)*h1 + eigv(i,2,2)*h2
531
532C Deviatoric part
533 fac = third*(hp1*ec(i,1) + hp2*ec(i,2) + h3*ec(i,3))
534 hd1 = hp1 - fac/max(em20,ec(i,1))
535 hd2 = hp2 - fac/max(em20,ec(i,2))
536 hd3 = h3 - fac/max(em20,ec(i,3))
537C add the viscoous stress
538 s(i,1) = s(i,1) + gama(k)*rvd(i)*hd1
539 s(i,2) = s(i,2) + gama(k)*rvd(i)*hd2
540 s(i,3) = s(i,3) + gama(k)*rvd(i)*hd3
541 ENDDO ! 1,NEL
542 ENDDO ! NPRONY
543
544C transformation PK2 to cauchy stress
545 DO i=1,nel
546 invr(i)= one/max(em20,rv(i))
547 t(i,1) = invr(i)*s(i,1)*ev(i,1)**2
548 t(i,2) = invr(i)*s(i,2)*ev(i,2)**2
549 t(i,3) = invr(i)*s(i,3)*ev(i,3)**2
550 ENDDO
551 ENDIF
552C-----------------------------------------------------------
553C-----------------------------------------------------------
554C Principal Cauchy viscous stress -> global directions
555 DO i=1,nel
556 signxx(i) = eigv(i,1,1)*t(i,1) + eigv(i,1,2)*t(i,2)
557 signyy(i) = eigv(i,2,1)*t(i,1) + eigv(i,2,2)*t(i,2)
558 signxy(i) = eigv(i,3,1)*t(i,1) + eigv(i,3,2)*t(i,2)
559 signyz(i) = sigoyz(i) + gs(i)*depsyz(i)
560 signzx(i) = sigozx(i) + gs(i)*depszx(i)
561 ENDDO
562C-----------------------------------------------------------
563C set sound speed & viscosity
564 DO i=1,nel
565 dezz(i) =-nu/(one-nu)*(depsxx(i)+depsyy(i))
566 thkn(i) = thkn(i) + dezz(i)*thklyl(i)
567 rho(i) = rho0(i)/rv(i)
568C
569 emax = gmax*(one + nu)
570 a11 = emax/(one - nu**2)
571 soundsp(i)= sqrt(a11/rho0(i))
572!! SOUNDSP(I) = SQRT((TWO_THIRD*GMAX+RBULK)/RHO(I))
573 viscmax(i) = zero
574 ENDDO
575
576C-----------
577 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
max
#define max(a, b)
Definition macros.h:21