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

Go to the source code of this file.

Functions/Subroutines

subroutine plas24b (nel, nindx, indx, ngl, pm, sig, dam, crak, epsvp, cdam, rho, eint, vk0, vk, rob, pla, deps1, deps2, deps3, deps4, deps5, deps6, s1, s2, s3, s4, s5, s6, scal1, scal2, scal3, scle2)

Function/Subroutine Documentation

◆ plas24b()

subroutine plas24b ( integer nel,
integer nindx,
integer, dimension(nindx) indx,
integer, dimension(nel) ngl,
pm,
dimension(nel,6) sig,
dimension(nel,3) dam,
dimension(nel,3) crak,
dimension(nel) epsvp,
dimension(nel,3,3) cdam,
dimension(nel) rho,
dimension(nel) eint,
dimension(nel) vk0,
dimension(nel) vk,
dimension(nel) rob,
dimension(nel,7) pla,
dimension(nel) deps1,
dimension(nel) deps2,
dimension(nel) deps3,
dimension(nel) deps4,
dimension(nel) deps5,
dimension(nel) deps6,
dimension(nel) s1,
dimension(nel) s2,
dimension(nel) s3,
dimension(nel) s4,
dimension(nel) s5,
dimension(nel) s6,
dimension(nel) scal1,
dimension(nel) scal2,
dimension(nel) scal3,
dimension(nel) scle2 )

Definition at line 28 of file plas24b.F.

34C-----------------------------------------------
35C I m p l i c i t T y p e s
36C-----------------------------------------------
37#include "implicit_f.inc"
38C-----------------------------------------------
39C C o m m o n B l o c k s
40C-----------------------------------------------
41#include "param_c.inc"
42#include "comlock.inc"
43C-----------------------------------------------
44C D u m m y A r g u m e n t s
45C-----------------------------------------------
46 INTEGER NEL,NINDX,NGL(NEL),INDX(NINDX)
47 my_real pm(npropm)
48 my_real, DIMENSION(NEL) :: s1,s2,s3,s4,s5,s6,epsvp,rho,eint,vk0,vk,rob,
49 . scal1,scal2,scal3,scle2,deps1,deps2,deps3,deps4,deps5,deps6
50 my_real, DIMENSION(NEL,7) :: pla
51 my_real, DIMENSION(NEL,3,3) :: cdam
52 my_real, DIMENSION(NEL,6) :: sig
53 my_real, DIMENSION(NEL,3) :: dam,crak
54C-----------------------------------------------
55C L o c a l V a r i a b l e s
56C-----------------------------------------------
57 INTEGER I,J,NITER,ITER
58 my_real, DIMENSION(NEL) :: sm,de1,de2,de3,de4,de5,de6,c44,c55,c66
59 my_real fc, rt, rc, rct1, rct2, aa, bc,
60 . bt, ac, hbp, ali0, alf0, vky, tol, rok0, ro0, hv0, vmax, expo,
61 . young, nu, g, den, rate, fac, fac1,rok, r2, aj3, cs3t,
62 . bb, df, rf, rf2, aj2, ajj, dkdsm, drfdsm, b0, drf3, b1,
63 . b2, alpha, phi, to, dfdto, hpv, hp, ecr,sigm,
64 . ts1, ts2, ts3, ts4, ts5, ts6,
65 . dfs1,dfs2,dfs3,dfs4,dfs5,dfs6,dgs1,dgs2,dgs3,dgs4,dgs5,dgs6,
66 . h1a, h2a, h3a, h4a,h5a, h6a, h1n, h2n, h3n, h4n, h5n, h6n, hh,
67 . lambda,depsp1,depsp2,depsp3,depsp4,depsp5,depsp6,
68 . depsel1,depsel2,depsel3,depsel4,depsel5,depsel6,
69 . depsvp, rr,vkmax, ro, div, vkk, numer, denom,
70 . bulk,dp,rho0,hvfac,d1,d2,d3,depsv,seuil,dfdj,norms,
71 . cp11,cp12,cp13,cp14,cp15,cp16,cp21,cp22,cp23,cp24,cp25,cp26,cp31,
72 . cp32,cp33,cp34,cp35,cp36,cp41,cp42,cp43,cp44,cp45,cp46,cp51,cp52,
73 . cp53,cp54,cp55,cp56,cp61,cp62,cp63,cp64,cp65,cp66,
74 . c11,c12,c13,c22,c23,c33
75C=======================================================================
76 lambda = zero
77 rho0 = pm(1)
78 young = pm(20)
79 nu = pm(21)
80 g = pm(22)
81 bulk = pm(32)
82 fc = pm(33)
83 rt = pm(34)
84 rc = pm(35)
85 rct1 = pm(36)
86 rct2 = pm(37)
87 aa = pm(38)
88 bc = pm(39)
89 bt = pm(40)
90 ac = pm(41)
91 hbp = pm(43)
92 ali0 = pm(44)
93 alf0 = pm(45)
94 vky = pm(46)
95 rok0 = pm(29)
96 ro0 = pm(30)
97 hv0 = pm(48)
98 vmax = pm(27)
99 expo = pm(49)
100 hvfac = pm(58)
101c
102 tol =( rt-rc)/twenty
103c------------------------------------------------
104C RETOUR SUR LE CRITERE ,SIG ,CDAM,CRAK,...
105c------------------------------------------------
106 DO j = 1,nindx
107 i = indx(j)
108 crak(i,1)=crak(i,1)-scle2(i)*deps1(i)
109 crak(i,2)=crak(i,2)-scle2(i)*deps2(i)
110 crak(i,3)=crak(i,3)-scle2(i)*deps3(i)
111 de1(i)=one - max( zero , sign(dam(i,1),crak(i,1)) )
112 de2(i)=one - max( zero , sign(dam(i,2),crak(i,2)) )
113 de3(i)=one - max( zero , sign(dam(i,3),crak(i,3)) )
114 scal1(i)=half+sign(half,de1(i)-one)
115 scal2(i)=half+sign(half,de2(i)-one)
116 scal3(i)=half+sign(half,de3(i)-one)
117 de4(i)=scal1(i)*scal2(i)
118 de5(i)=scal2(i)*scal3(i)
119 de6(i)=scal3(i)*scal1(i)
120c-------
121 den = one - nu**2 *(de4(i) + de5(i) + de6(i)
122 . + two*nu*scal1(i)*scal2(i)*scal3(i))
123 den = one / den
124 cdam(i,1,1) = young*de1(i)*(one -nu**2*scal2(i)*scal3(i))*den
125 cdam(i,2,2) = young*de2(i)*(one -nu**2*scal1(i)*scal3(i))*den
126 cdam(i,3,3) = young*de3(i)*(one -nu**2*scal2(i)*scal1(i))*den
127 cdam(i,1,2) = nu*young*scal1(i)*scal2(i) *(one+nu*scal3(i))*den
128 cdam(i,2,3) = nu*young*scal2(i)*scal3(i) *(one+nu*scal1(i))*den
129 cdam(i,1,3) = nu*young*scal1(i)*scal3(i) *(one+nu*scal2(i))*den
130 cdam(i,2,1) = cdam(i,1,2)
131 cdam(i,3,1) = cdam(i,1,3)
132 cdam(i,3,2) = cdam(i,2,3)
133 c44(i) = g*de4(i)
134 c55(i) = g*de5(i)
135 c66(i) = g*de6(i)
136c-----
137 sig(i,1)=cdam(i,1,1)*crak(i,1)+cdam(i,1,2)*crak(i,2)+cdam(i,1,3)*crak(i,3)
138 sig(i,2)=cdam(i,2,1)*crak(i,1)+cdam(i,2,2)*crak(i,2)+cdam(i,2,3)*crak(i,3)
139 sig(i,3)=cdam(i,3,1)*crak(i,1)+cdam(i,3,2)*crak(i,2)+cdam(i,3,3)*crak(i,3)
140 sig(i,4)=de4(i)*sig(i,4)-c44(i)*deps4(i)*(one-scle2(i))
141 sig(i,5)=de5(i)*sig(i,5)-c55(i)*deps5(i)*(one-scle2(i))
142 sig(i,6)=de6(i)*sig(i,6)-c66(i)*deps6(i)*(one-scle2(i))
143 s1(i) = scal1(i)*sig(i,1)
144 s2(i) = scal2(i)*sig(i,2)
145 s3(i) = scal3(i)*sig(i,3)
146 s4(i) = sig(i,4)*scal1(i)*scal2(i)
147 s5(i) = sig(i,5)*scal2(i)*scal3(i)
148 s6(i) = sig(i,6)*scal3(i)*scal1(i)
149 sm(i) = (s1(i)+s2(i)+s3(i)) * third
150 s1(i) = s1(i) - sm(i)
151 s2(i) = s2(i) - sm(i)
152 s3(i) = s3(i) - sm(i)
153c----------------------------
154 denom = abs(crak(i,1))+abs(crak(i,2))+abs(crak(i,3))
155 . +(abs(sig(i,4)) +abs(sig(i,5)) +abs(sig(i,6)))/g
156c----------------------------
157 IF (denom == zero) cycle
158c----------------------------
159 numer = abs(deps1(i))+abs(deps2(i))+abs(deps3(i))
160 . + abs(deps4(i))+abs(deps5(i))+abs(deps6(i))
161c
162 rate = numer / denom
163 niter = nint(three*rate) + 1
164 niter = min0(niter,10)
165c--------------------------------
166 fac = scle2(i)/niter
167 deps1(i) = fac * deps1(i)
168 deps2(i) = fac * deps2(i)
169 deps3(i) = fac * deps3(i)
170 deps4(i) = fac * deps4(i)
171 deps5(i) = fac * deps5(i)
172 deps6(i) = fac * deps6(i)
173c---------------------------------------------------------
174c ... DEBUT DES ITERATIONS
175c---------------------------------------------------------
176 DO iter=1,niter
177c--------------------
178 rok = rok0+rob(i)-ro0
179 r2 = s1(i)**2 + s2(i)**2 + s3(i)**2
180 . +(s4(i)**2 + s5(i)**2 + s6(i)**2)*two
181 IF (sm(i) >= rt-tol) THEN
182 vk(i) = one
183 ELSEIF (sm(i) > rc) THEN
184 vk(i) = one +(one-vk0(i))*(rct1-two*rc*sm(i)+sm(i)**2)/rct2
185 ELSEIF (sm(i) > rok)THEN
186 vk(i) = vk0(i)
187C CAPE
188 ELSEIF (sm(i) > rob(i)) THEN
189 vk(i) = vk0(i)*(one - ((sm(i)-rok)/(rob(i)-rok))**2)
190 ELSE
191 vk(i) = zero
192 ENDIF
193C ON TRONQUE PRES DE l'AXE des pressions
194C CRITERE
195 aj3 = s1(i)*s2(i)*s3(i)
196 . - s1(i)*s5(i)*s5(i)-s2(i)*s6(i)*s6(i)-s3(i)*s4(i)*s4(i)
197 . + two*s4(i)*s5(i)*s6(i)
198 cs3t= half * aj3*(three/half*max(r2,em20))**three_half
199 cs3t= min(one,cs3t)
200 cs3t= max(-one,cs3t)
201 bb = half*((one - cs3t)*bc+(one +cs3t)*bt)
202 df = sqrt(bb*bb + max(-aa*sm(i) + ac,em9))
203 rf = (-bb+df)/aa
204 rf2 = rf**2
205 aj2 = half*r2
206 ajj = sqrt(aj2)
207 rr = sqr2*ajj
208c
209 IF (r2/rf2 > em04) THEN
210C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
211C NOW COMPUTES ALL TERMS TO BUILD UP PLASTIC MATRIX
212C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
213C D I M E N S I O N S RF:CONT TO:CONT PHI:ADIM B0:ADIM B1:1/CONT
214C B2:1/CONT**2 DFDTO:ADIM DRFDSM:ADIM DKDSM:1/CONT
215C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
216C PRINCIPALES DERIVEES
217 sm(i) = max(rob(i),sm(i))
218 IF (sm(i) >= rt-tol) THEN
219 dkdsm =zero
220 ELSEIF(sm(i)>rc)THEN
221 dkdsm = two*(one -vk0(i))*(sm(i)-rc)/rct2
222 ELSEIF(sm(i)>rok)THEN
223 dkdsm = zero
224 ELSE
225 dkdsm = -two*vk0(i)*(sm(i)-rok)/(rob(i)-rok)**2
226 ENDIF
227 drfdsm= -half / df
228 b0 = -third * vk(i) *drfdsm - third * rf * dkdsm
229 drf3 = half* (-one + bb/df) * (bt-bc)/aa
230 b1 = half*sqr2/ajj
231 . + vk(i)*drf3*fourth*aj3*(three/aj2)**twop5
232 b2 = -vk(i)*drf3*half*(three/aj2)**three_half
233C DF/DSIG
234 ts1 = s1(i)**2 + s4(i)**2 + s6(i)**2 - two*third*aj2
235 ts2 = s2(i)**2 + s4(i)**2 + s5(i)**2 - two*third*aj2
236 ts3 = s3(i)**2 + s5(i)**2 + s6(i)**2 - two*third*aj2
237 ts4 = two* (s5(i)*s6(i)-s4(i)*s3(i))
238 ts5 = two* (s6(i)*s4(i)-s5(i)*s1(i))
239 ts6 = two* (s4(i)*s5(i)-s6(i)*s2(i))
240 dfs1= b0 + b1*s1(i) + b2*ts1
241 dfs2= b0 + b1*s2(i) + b2*ts2
242 dfs3= b0 + b1*s3(i) + b2*ts3
243 dfs4= two*b1*s4(i) + b2*ts4
244 dfs5= two*b1*s5(i) + b2*ts5
245 dfs6= two*b1*s6(i) + b2*ts6
246C
247 IF (sm(i) > rok .and. vk(i) > vky) THEN
248 alpha = (one-vk(i))*ali0 + (vk(i)-vky)*alf0
249 alpha = alpha / (one-vky)
250 ELSEIF (vk0(i) > vky) THEN
251 alpha = (one-vk0(i))*ali0 + (vk0(i)-vky)*alf0
252 alpha = alpha / (one-vky)
253 ELSE
254 alpha = ali0
255 ENDIF
256C
257c Definition de la zone de transition
258C
259c SEUIL = ALI0
260 seuil = min(-em01,ali0)
261c IF(SM(I)<ROK )THEN
262c FAC=(SEUIL-B0)/(SEUIL+THIRD * VK0(I) *DRFDSM)
263 IF (sm(i) < rok .AND. b0 < zero) THEN
264 fac = (seuil-b0) / seuil
265 fac = max(zero,fac)
266 fac = min(one,fac)
267 fac1= one-fac
268 ELSE
269 fac = one
270 fac1= zero
271 ENDIF
272c Securite on ne veut pas que la normale a la loi d'ecoulement soit rentrante
273c + on liminte le cos de l'angle max entre citere etloi d'ecoulement a 0.5
274c dans le repere I1 J=sqrt(J2) ng=normale a G (alpha,1)
275c normal a f nf=(B0,DFDJ)
276c cos= ng.nf/||nf||||ng||
277 dfdj = two*b1*ajj
278cc NORMS=SQRT((ALPHA**2+ONE)*(B0**2+DFDJ**2))
279 norms=sqrt((b0**2+dfdj**2))
280 IF (b0 >=zero) THEN
281 alpha = max(alpha, (em01*norms - dfdj)/max(em20,b0))
282 ELSE
283 alpha = min(alpha, (em01*norms - dfdj)/min(b0,-em20))
284 ENDIF
285c Dilatation max atteinte
286 IF (rho(i) < rho0 .OR. eint(i)<0) alpha=zero
287c Ecoulement non associe
288 dgs1=fac*alpha+s1(i)/(two*ajj)
289 dgs2=fac*alpha+s2(i)/(two*ajj)
290 dgs3=fac*alpha+s3(i)/(two*ajj)
291 dgs4=s4(i)/ajj
292 dgs5=s5(i)/ajj
293 dgs6=s6(i)/ajj
294c Transition et associe
295 dgs1=fac*dgs1+fac1*dfs1
296 dgs2=fac*dgs2+fac1*dfs2
297 dgs3=fac*dgs3+fac1*dfs3
298 dgs4=fac*dgs4+fac1*dfs4
299 dgs5=fac*dgs5+fac1*dfs5
300 dgs6=fac*dgs6+fac1*dfs6
301c compactance dilatance effective (!= alpha sur la cape)
302c correction erreur ALPHA = DGS1+DGS2+DGS3
303 alpha = third*(dgs1+dgs2+dgs3)
304C ECROUISSAGE
305c contribution deviatorique
306 IF (sm(i) < rok) THEN
307 to = sqr3_2*vk0(i)*rf
308 ELSE
309 to = sqr3_2*vk(i)*rf
310 ENDIF
311 phi = max(zero,alpha*three*sm(i)+ajj)/to
312 dfdto=-sqrt(two_third)+b0
313 ecr = phi*dfdto*hbp*fac
314c ECR= TEN*PHI*DFDTO*HBP*FAC*(ONE-VK0(I))/(ONE-VKY) ! ecrousissage progressif
315c IF (VK0(I)==ONE ) ECR=ZERO
316c Contribution volumetrique
317 hpv =hv0*(one-hvfac*vk(i))*max(one,exp(expo*epsvp(i)))
318 IF (sm(i) >= rok) THEN
319 IF(alpha>=zero) hpv=zero
320 ELSE
321 IF (alpha>=zero) THEN ! si dilatant on ne bouge pas la cape
322 hpv = zero
323 ELSE
324 ecr = ecr + rf * dkdsm*hpv*alpha
325 ENDIF
326 ENDIF
327c plasticite parfaite sur le critere de rupture
328c ECR=ECR*(HALF-SIGN(HALF,VK(I)-ONE))
329c write(993,"(A,3E12.3)")'vkvk0 ',VK(I),VK0(I),SQRT(R2/RF2)
330c write(993,"(A,3E12.3)")'SMrobrk',sm(I),ROB(I),ROK
331 IF(r2>=rf2 .OR.vk(i)>=one )THEN
332 hpv=zero
333 ecr=zero
334 ENDIF
335C PRODUITS TENSORIELS
336C *scal car si la fissure est ouverte pas de plasticite dans cette direction
337C les termes croises sont deja nul dans cdam
338 dfs1=dfs1*scal1(i)
339 dfs2=dfs2*scal2(i)
340 dfs3=dfs3*scal3(i)
341 dfs4=dfs4*de4(i)
342 dfs5=dfs5*de5(i)
343 dfs6=dfs6*de6(i)
344 dgs1=dgs1*scal1(i)
345 dgs2=dgs2*scal2(i)
346 dgs3=dgs3*scal3(i)
347 dgs4=dgs4*de4(i)
348 dgs5=dgs5*de5(i)
349 dgs6=dgs6*de6(i)
350C
351 h1a = cdam(i,1,1)*dfs1 + cdam(i,2,1)*dfs2 + cdam(i,3,1)*dfs3
352 h2a = cdam(i,1,2)*dfs1 + cdam(i,2,2)*dfs2 + cdam(i,3,2)*dfs3
353 h3a = cdam(i,1,3)*dfs1 + cdam(i,2,3)*dfs2 + cdam(i,3,3)*dfs3
354 h4a = c44(i)*dfs4
355 h5a = c55(i)*dfs5
356 h6a = c66(i)*dfs6
357C
358 h1n = cdam(i,1,1)*dgs1 + cdam(i,1,2)*dgs2 + cdam(i,1,3)*dgs3
359 h2n = cdam(i,2,1)*dgs1 + cdam(i,2,2)*dgs2 + cdam(i,2,3)*dgs3
360 h3n = cdam(i,3,1)*dgs1 + cdam(i,3,2)*dgs2 + cdam(i,3,3)*dgs3
361 h4n = c44(i)*dgs4
362 h5n = c55(i)*dgs5
363 h6n = c66(i)*dgs6
364C LAMBDA
365 hh= dfs1*h1n+dfs2*h2n+dfs3*h3n+dfs4*h4n+dfs5*h5n+dfs6*h6n
366 . - min(zero,ecr)
367C
368 lambda= h1a*deps1(i)+h2a*deps2(i)+h3a*deps3(i)
369 . +h4a*deps4(i)+h5a*deps5(i)+h6a*deps6(i)
370 lambda=lambda/hh
371!---
372! avoid negative plastic strain
373 lambda=max(lambda,zero)
374!---
375C DEFORMATIONS PLASTIQUES
376 depsp1=lambda*dgs1
377 depsp2=lambda*dgs2
378 depsp3=lambda*dgs3
379 depsp4=lambda*dgs4
380 depsp5=lambda*dgs5
381 depsp6=lambda*dgs6
382 depsvp=depsp1+depsp2+depsp3
383 epsvp(i)=epsvp(i)+min(zero,depsvp) ! on cumule seulement la compaction
384c
385C DEFORMATIONS ELASTIQUES
386 depsel1=deps1(i)-depsp1
387 depsel2=deps2(i)-depsp2
388 depsel3=deps3(i)-depsp3
389 depsel4=deps4(i)-depsp4
390 depsel5=deps5(i)-depsp5
391 depsel6=deps6(i)-depsp6
392
393c SIG(I,1)=SIG(I,1)+CDAM(I,1,1)*DEPSEL1+CDAM(I,1,2)*DEPSEL2
394c . +CDAM(I,1,3)*DEPSEL3
395c SIG(I,2)=SIG(I,2)+CDAM(I,2,1)*DEPSEL1+CDAM(I,2,2)*DEPSEL2
396c . +CDAM(I,2,3)*DEPSEL3
397c SIG(I,3)=SIG(I,3)+CDAM(I,3,1)*DEPSEL1+CDAM(I,3,2)*DEPSEL2
398c . +CDAM(I,3,3)*DEPSEL3
399c SIG(I,4)=SIG(I,4)+C44(I)*DEPSEL4
400c SIG(I,5)=SIG(I,5)+C55(I)*DEPSEL5
401c SIG(I,6)=SIG(I,6)+C66(I)*DEPSEL6
402C
403C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
404C RECALCUL DES DEFORMATION ELASTIQUES POUR REOUVERTURE DES FISSURES
405C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
406c C11 = ONE/DE1(I)/YOUNG
407c C12 = -NU*SCAL1(I)*SCAL2(I)/YOUNG
408c C13 = -NU*SCAL1(I)*SCAL3(I)/YOUNG
409c CRAK(I,1) = C11*SIG(I,1)+C12*SIG(I,2)+C13*SIG(I,3)
410c C22 = ONE/DE2(I)/YOUNG
411c C23 = -NU*SCAL2(I)*SCAL3(I)/YOUNG
412c CRAK(I,2) = C12*SIG(I,1)+C22*SIG(I,2)+C23*SIG(I,3)
413c C33 = ONE/DE3(I)/YOUNG
414c CRAK(I,3) = C13*SIG(I,1)+C23*SIG(I,2)+C33*SIG(I,3)
415C
416 crak(i,1) =crak(i,1)+depsel1
417 crak(i,2) =crak(i,2)+depsel2
418 crak(i,3) =crak(i,3)+depsel3
419c
420 de1(i)=one - max( zero , sign(dam(i,1),crak(i,1)) )
421 de2(i)=one - max( zero , sign(dam(i,2),crak(i,2)) )
422 de3(i)=one - max( zero , sign(dam(i,3),crak(i,3)) )
423 scal1(i)=half+sign(half,de1(i)-one)
424 scal2(i)=half+sign(half,de2(i)-one)
425 scal3(i)=half+sign(half,de3(i)-one)
426 de4(i)=scal1(i)*scal2(i)
427 de5(i)=scal2(i)*scal3(i)
428 de6(i)=scal3(i)*scal1(i)
429c--------------------------
430 den = one - nu**2 *(de4(i) + de5(i) + de6(i)
431 . + two*nu*scal1(i)*scal2(i)*scal3(i))
432C
433 cdam(i,1,1) = young*de1(i)*(one-nu**2*scal2(i)*scal3(i))/den
434 cdam(i,2,2) = young*de2(i)*(one -nu**2*scal1(i)*scal3(i))/den
435 cdam(i,3,3) = young*de3(i)*(one -nu**2*scal2(i)*scal1(i))/den
436 cdam(i,1,2) = nu*young*scal1(i)*scal2(i) *(one+nu*scal3(i))/den
437 cdam(i,2,3) = nu*young*scal2(i)*scal3(i) *(one+nu*scal1(i))/den
438 cdam(i,1,3) = nu*young*scal1(i)*scal3(i) *(one+nu*scal2(i))/den
439 cdam(i,2,1) = cdam(i,1,2)
440 cdam(i,3,1) = cdam(i,1,3)
441 cdam(i,3,2) = cdam(i,2,3)
442 sig(i,1) = cdam(i,1,1)*crak(i,1)+cdam(i,1,2)*crak(i,2)+cdam(i,1,3)*crak(i,3)
443 sig(i,2) = cdam(i,2,1)*crak(i,1)+cdam(i,2,2)*crak(i,2)+cdam(i,2,3)*crak(i,3)
444 sig(i,3) = cdam(i,3,1)*crak(i,1)+cdam(i,3,2)*crak(i,2)+cdam(i,3,3)*crak(i,3)
445 sig(i,4) = (sig(i,4)+c44(i)*depsel4)*de4(i)
446 sig(i,5) = (sig(i,5)+c55(i)*depsel5)*de5(i)
447 sig(i,6) = (sig(i,6)+c66(i)*depsel6)*de6(i)
448c SIG(I,4) = SIG(I,4)*DE4(I)
449c SIG(I,5) = SIG(I,5)*DE5(I)
450c SIG(I,6) = SIG(I,6)*DE6(I)
451c------------------------------------------------
452 s1(i) = scal1(i)*sig(i,1)
453 s2(i) = scal2(i)*sig(i,2)
454 s3(i) = scal3(i)*sig(i,3)
455 s4(i) = sig(i,4)*scal1(i)*scal2(i)
456 s5(i) = sig(i,5)*scal2(i)*scal3(i)
457 s6(i) = sig(i,6)*scal3(i)*scal1(i)
458 sm(i) = third * (s1(i)+s2(i)+s3(i))
459C
460 s1(i) = s1(i)-sm(i)
461 s2(i) = s2(i)-sm(i)
462 s3(i) = s3(i)-sm(i)
463C . . . . . . . . . . . . . . . . . . . . . . . . . . . .
464C CALCUL DES PARAMETRES D' ECROUISSAGE COURANTS VK0,ROB
465C CORRECTIONS EVENTUELLES
466C . . . . . . . . . . . . . . . . . . . . . . . . . . . .
467 IF (sm(i) > ac/aa)
468 . sm(i)=sm(i)-three*(sm(i)-ac/aa)/(scal1(i)+scal2(i)+scal3(i))
469C
470 r2 = s1(i)**2 + s2(i)**2 + s3(i)**2
471 . +(s4(i)**2 + s5(i)**2 + s6(i)**2)*two
472 rr = sqrt(r2)
473 aj3 = s1(i)*s2(i)*s3(i)
474 . - s1(i)*s5(i)*s5(i) - s2(i)*s6(i)*s6(i) - s3(i)*s4(i)*s4(i)
475 . + two*s4(i)*s5(i)*s6(i)
476 cs3t= half * aj3*(three/half*max(r2,em20))**three_half
477 cs3t= min(one,cs3t)
478 cs3t= max(-one,cs3t)
479 bb = half*((one-cs3t)*bc + (one+cs3t)*bt)
480 df = sqrt(bb*bb + max(-aa*sm(i)+ac, zero))
481 rf = (-bb + df) / aa
482c
483 IF (rf == zero) THEN
484 vk(i) = zero
485 ELSE
486 vk(i) = rr / rf
487 ENDIF
488C
489c on reprojette si en dehors du critere
490 IF (vk(i) > one) THEN
491 fac = one/vk(i)
492 IF (scal1(i) > zep9) sig(i,1) = s1(i)*fac+sm(i)
493 IF (scal2(i) > zep9) sig(i,2) = s2(i)*fac+sm(i)
494 IF (scal3(i) > zep9) sig(i,3) = s3(i)*fac+sm(i)
495 IF (scal1(i)*scal2(i) > zep9) sig(i,4) = s4(i)*fac
496 IF (scal2(i)*scal3(i) > zep9) sig(i,5) = s5(i)*fac
497 IF (scal3(i)*scal1(i) > zep9) sig(i,6) = s6(i)*fac
498 vk(i) = one
499 c11 = one/de1(i)/young
500 c12 = -nu*scal1(i)*scal2(i)/young
501 c13 = -nu*scal1(i)*scal3(i)/young
502 crak(i,1) = c11*sig(i,1)+c12*sig(i,2)+c13*sig(i,3)
503 c22 = one/de2(i)/young
504 c23 = -nu*scal2(i)*scal3(i)/young
505 crak(i,2) = c12*sig(i,1)+c22*sig(i,2)+c23*sig(i,3)
506 c33 = one/de3(i)/young
507 crak(i,3) = c13*sig(i,1)+c23*sig(i,2)+c33*sig(i,3)
508C
509 de1(i)=one - max( zero , sign(dam(i,1),crak(i,1)) )
510 de2(i)=one - max( zero , sign(dam(i,2),crak(i,2)) )
511 de3(i)=one - max( zero , sign(dam(i,3),crak(i,3)) )
512 scal1(i)=half + sign(half,de1(i)-one)
513 scal2(i)=half + sign(half,de2(i)-one)
514 scal3(i)=half + sign(half,de3(i)-one)
515 de4(i)=scal1(i)*scal2(i)
516 de5(i)=scal2(i)*scal3(i)
517 de6(i)=scal3(i)*scal1(i)
518c---------------------------------
519 den = one - nu**2 *(de4(i) + de5(i) + de6(i)
520 . + two*nu*scal1(i)*scal2(i)*scal3(i))
521C
522 cdam(i,1,1) = young*de1(i)*(one -nu**2*scal2(i)*scal3(i))/den
523 cdam(i,2,2) = young*de2(i)*(one -nu**2*scal1(i)*scal3(i))/den
524 cdam(i,3,3) = young*de3(i)*(one -nu**2*scal2(i)*scal1(i))/den
525 cdam(i,1,2) = nu*young*scal1(i)*scal2(i) *(one +nu*scal3(i))/den
526 cdam(i,2,3) = nu*young*scal2(i)*scal3(i) *(one +nu*scal1(i))/den
527 cdam(i,1,3) = nu*young*scal1(i)*scal3(i) *(one +nu*scal2(i))/den
528 cdam(i,2,1) = cdam(i,1,2)
529 cdam(i,3,1) = cdam(i,1,3)
530 cdam(i,3,2) = cdam(i,2,3)
531 sig(i,1)=cdam(i,1,1)*crak(i,1)+cdam(i,1,2)*crak(i,2)+cdam(i,1,3)*crak(i,3)
532 sig(i,2)=cdam(i,2,1)*crak(i,1)+cdam(i,2,2)*crak(i,2)+cdam(i,2,3)*crak(i,3)
533 sig(i,3)=cdam(i,3,1)*crak(i,1)+cdam(i,3,2)*crak(i,2)+cdam(i,3,3)*crak(i,3)
534 sig(i,4)=de4(i)*sig(i,4)
535 sig(i,5)=de5(i)*sig(i,5)
536 sig(i,6)=de6(i)*sig(i,6)
537 ENDIF
538C
539 rob(i) = min(rob(i), rob(i)+hpv*depsvp, sm(i))
540 rok = rob(i) - ro0 + rok0
541
542 IF (sm(i) >= rt-tol)THEN
543 vkk=one
544 ELSEIF(sm(i) > rc) THEN
545 div= min(-em10,rct1- two*rc*sm(i)+sm(i)*sm(i) )
546 vkk=one +(one - vk(i))*rct2/div
547 ELSEIF(sm(i)>=rok) THEN
548 vkk=vk(i)
549 ELSE
550c si au dela de la cape, on change rob
551 vkmax=one- ((sm(i)-rok)/(rob(i)-rok))**2
552 IF (vk(i)>vkmax) THEN
553 rok = sm(i) - sqrt(one-vk(i))*(ro0-rok0)
554 rob(i) = rok+ro0-rok0
555 vkk = one
556 ELSE
557c sinon on calcule k0
558 vkk = vk(i)/max(vkmax,em03)
559 ENDIF
560 ENDIF
561 vkk = min(vkk,one)
562 vk0(i)= max(vkk,vk0(i))
563c write(7,*)'normal'
564c
565 ELSE
566c PURE TRIAXIAL TO SOLVE UNDETERMINATION
567 depsv = (deps1(i) + deps2(i) + deps3(i))
568 d1 = deps1(i) - depsv*third
569 d2 = deps2(i) - depsv*third
570 d3 = deps3(i) - depsv*third
571 hpv = hv0*max(one, exp(expo*epsvp(i)))
572 depsvp = bulk/(bulk+hpv)*depsv
573 epsvp(i)= epsvp(i) + depsvp
574 ro = min( rob(i),rob(i) + hpv*depsvp )
575 dp = ro - rob(i)
576 sig(i,1) = sig(i,1) + dp + two*g*d1
577 sig(i,2) = sig(i,2) + dp + two*g*d2
578 sig(i,3) = sig(i,3) + dp + two*g*d3
579 sig(i,4) = sig(i,4) + g*deps4(i)
580 sig(i,5) = sig(i,5) + g*deps5(i)
581 sig(i,6) = sig(i,6) + g*deps6(i)
582 rob(i) = ro
583 vk0(i) = one
584c write(7,*)'pure triaxial'
585 ENDIF
586c equivalent plastic strain for output
587 r2 = (s1(i)**2 + s2(i)**2 + s3(i)**2
588 . + (s4(i)**2 + s5(i)**2 + s6(i)**2)*two)*three_half
589 sigm = sig(i,1) + sig(i,2) + sig(i,2)
590 rr = one/sqr3
591 IF (r2 > zero) rr = rr + alpha*sigm/sqrt(r2)
592!! PLA(I) = PLA(I) + LAMBDA * RR
593 pla(i,1) = pla(i,1) + lambda * rr ! --- plastic strain --- scalar
594!---
595! full plastic strain tensor (for output only)
596!---
597 pla(i,2) = pla(i,2) + lambda*dgs1
598 pla(i,3) = pla(i,3) + lambda*dgs2
599 pla(i,4) = pla(i,4) + lambda*dgs3
600 pla(i,5) = pla(i,5) + lambda*dgs4
601 pla(i,6) = pla(i,6) + lambda*dgs5
602 pla(i,7) = pla(i,7) + lambda*dgs6
603c-----------
604c FIN DES ITERATIONS
605 ENDDO ! ITER
606c-----------
607 ENDDO ! NINDX
608c-----------
609 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
alpha
#define alpha
Definition eval.h:35
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21