OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
plas24.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!|| plas24 ../engine/source/materials/mat/mat024/plas24.F
25!||--- called by ------------------------------------------------------
26!|| conc24 ../engine/source/materials/mat/mat024/conc24.F
27!||====================================================================
28 SUBROUTINE plas24(NEL ,NINDX ,INDX ,NGL ,PM ,
29 . SIG ,DAM ,CRAK ,
30 . RHO ,EINT ,VK0 ,VK ,ROB ,CDAM ,
31 . E1 ,E2 ,E3 ,E4 ,E5 ,E6 ,
32 . S1 ,S2 ,S3 ,S4 ,S5 ,S6 ,
33 . SCAL1 ,SCAL2 ,SCAL3 ,SCLE2 )
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 NINDX,NEL,NGL(NEL),INDX(NINDX)
47 my_real PM(NPROPM)
48 my_real, DIMENSION(NEL) :: S1,S2,S3,S4,S5,S6,
49 . SCAL1,SCAL2,SCAL3,SCLE2,E1,E2,E3,E4,E5,E6,EINT,RHO,VK0,VK,ROB
50 my_real, DIMENSION(NEL,3,3) :: CDAM
51 my_real, DIMENSION(NEL,6) :: sig
52 my_real, DIMENSION(NEL,3) :: dam,crak
53C-----------------------------------------------
54C L o c a l V a r i a b l e s
55C-----------------------------------------------
56 INTEGER I,J,NITER,ITER,ICAP,IBUG
57 my_real, DIMENSION(NEL) :: sm,de1,de2,de3,de4,de5,de6,c44,c55,c66
58 my_real
59 . 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, rok, r2, aj3, cs3t,
62 . bb, df, rf, rf2, aj2, ajj, dkdsm, drfdsm, b0, drf3, b1,
63 . b2, alpha, phi, dfdto1, dfdto2, to, sq32, dfdto, hpv, hp, ecr,
64 . ts1, ts2, ts3, ts4, ts5, ts6, dfs1, dfs2, dfs3, dfs4, dfs5,
65 . dfs6, dgs1, dgs2, dgs3, dgs4, dgs5, dgs6, h1a, h2a, h3a, h4a,
66 . h5a, h6a, h1n, h2n, h3n, h4n, h5n, h6n, hh, cp11, cp12, cp13,
67 . cp14, cp15, cp16, cp21, cp22, cp23, cp24, cp25, cp26, cp31,
68 . cp32, cp33, cp34, cp35, cp36, cp41, cp42, cp43, cp44, cp45,
69 . cp46, cp51, cp52, cp53, cp54, cp55, cp56, cp61, cp62, cp63,
70 . cp64, cp65, cp66, c11, c12, c13, c22, c23, c33, vkold, rr,
71 . vkmax, ro, div, dvk, vkk, numer, denom,bulk,dp,rho0,dfdro
72C=======================================================================
73 rho0 = pm(1)
74 young = pm(20)
75 nu = pm(21)
76 g = pm(22)
77 vmax = pm(27)
78 rok0 = pm(29)
79 ro0 = pm(30)
80 bulk = pm(32)
81 fc = pm(33)
82 rt = pm(34)
83 rc = pm(35)
84 rct1 = pm(36)
85 rct2 = pm(37)
86 aa = pm(38)
87 bc = pm(39)
88 bt = pm(40)
89 ac = pm(41)
90 hbp = pm(43)
91 ali0 = pm(44)
92 alf0 = pm(45)
93 vky = pm(46)
94 hv0 = pm(48)
95 expo = pm(49)
96 icap = nint(pm(57))
97 ibug = nint(pm(59))
98c
99 tol = (rt - rc)/twenty
100c------------------------------------------------
101C RETOUR SUR LE CRITERE ,SIG ,CDAM,CRAK,...
102c------------------------------------------------
103 DO j = 1,nindx
104 i = indx(j)
105 crak(i,1) = crak(i,1)-scle2(i)*e1(i)
106 crak(i,2) = crak(i,2)-scle2(i)*e2(i)
107 crak(i,3) = crak(i,3)-scle2(i)*e3(i)
108 de1(i) = one - max( zero , sign(dam(i,1),crak(i,1)) )
109 de2(i) = one - max( zero , sign(dam(i,2),crak(i,2)) )
110 de3(i) = one - max( zero , sign(dam(i,3),crak(i,3)) )
111 scal1(i)= half+sign(half,de1(i)-one)
112 scal2(i)= half+sign(half,de2(i)-one)
113 scal3(i)= half+sign(half,de3(i)-one)
114 de4(i) = scal1(i)*scal2(i)
115 de5(i) = scal2(i)*scal3(i)
116 de6(i) = scal3(i)*scal1(i)
117c-------
118 den = one - nu**2 *(de4(i) + de5(i) + de6(i)
119 . + two*nu*scal1(i)*scal2(i)*scal3(i))
120 cdam(i,1,1) = young*de1(i)*(one -nu**2*scal2(i)*scal3(i))/den
121 cdam(i,2,2) = young*de2(i)*(one -nu**2*scal1(i)*scal3(i))/den
122 cdam(i,3,3) = young*de3(i)*(one -nu**2*scal2(i)*scal1(i))/den
123 cdam(i,1,2) = nu*young*scal1(i)*scal2(i) *(one+nu*scal3(i))/den
124 cdam(i,2,3) = nu*young*scal2(i)*scal3(i) *(one+nu*scal1(i))/den
125 cdam(i,1,3) = nu*young*scal1(i)*scal3(i) *(one+nu*scal2(i))/den
126 cdam(i,2,1) = cdam(i,1,2)
127 cdam(i,3,1) = cdam(i,1,3)
128 cdam(i,3,2) = cdam(i,2,3)
129 c44(i) = g*de4(i)
130 c55(i) = g*de5(i)
131 c66(i) = g*de6(i)
132c-----
133 sig(i,1)=cdam(i,1,1)*crak(i,1)+cdam(i,1,2)*crak(i,2)+cdam(i,1,3)*crak(i,3)
134 sig(i,2)=cdam(i,2,1)*crak(i,1)+cdam(i,2,2)*crak(i,2)+cdam(i,2,3)*crak(i,3)
135 sig(i,3)=cdam(i,3,1)*crak(i,1)+cdam(i,3,2)*crak(i,2)+cdam(i,3,3)*crak(i,3)
136 IF (ibug == 0) THEN
137 sig(i,4) = de4(i)*sig(i,4) - c44(i)*e4(i)*(one-scle2(i))
138 sig(i,5) = de5(i)*sig(i,5) - c55(i)*e5(i)*(one-scle2(i))
139 sig(i,6) = de6(i)*sig(i,6) - c66(i)*e6(i)*(one-scle2(i))
140 ELSE
141 sig(i,4) = de4(i)*sig(i,4) + c44(i)*e4(i)*scle2(i)
142 sig(i,5) = de5(i)*sig(i,5) + c55(i)*e5(i)*scle2(i)
143 sig(i,6) = de6(i)*sig(i,6) + c66(i)*e6(i)*scle2(i)
144 ENDIF
145 s1(i) = scal1(i)*sig(i,1)
146 s2(i) = scal2(i)*sig(i,2)
147 s3(i) = scal3(i)*sig(i,3)
148 s4(i) = sig(i,4)*scal1(i)*scal2(i)
149 s5(i) = sig(i,5)*scal2(i)*scal3(i)
150 s6(i) = sig(i,6)*scal3(i)*scal1(i)
151 sm(i) = third * (s1(i)+s2(i)+s3(i))
152 s1(i) = s1(i) - sm(i)
153 s2(i) = s2(i) - sm(i)
154 s3(i) = s3(i) - sm(i)
155C
156 numer = abs(e1(i))+abs(e2(i))+abs(e3(i))
157 . + abs(e4(i))+abs(e5(i))+abs(e6(i))
158c
159 denom = (abs(crak(i,1))+ abs(crak(i,2)) + abs(crak(i,3))
160 . + (abs(sig(i,4)) + abs(sig(i,5)) + abs(sig(i,6)))/g )
161c------------------------------------------------
162 IF (denom == zero) cycle
163c------------------------------------------------
164 rate = numer/denom
165 niter = nint(three*rate) + 1
166 niter = min0(niter,10)
167c------------------------------------------------
168 fac = scle2(i)/niter
169 e1(i) = fac * e1(i)
170 e2(i) = fac * e2(i)
171 e3(i) = fac * e3(i)
172 e4(i) = fac * e4(i)
173 e5(i) = fac * e5(i)
174 e6(i) = fac * e6(i)
175c---------------------------------------------------------
176c ... DEBUT DES ITERATIONS
177c---------------------------------------------------------
178 DO iter=1,niter
179c--------------------
180 rok = rok0+rob(i)-ro0
181 r2 = s1(i)**2+s2(i)**2+s3(i)**2
182 . + two*s4(i)**2+two*s5(i)**2+two*s6(i)**2
183 IF (sm(i) >= rt-tol)THEN
184 vk(i)=one
185 ELSEIF (sm(i) > rc) THEN
186 vk(i)=one +(one-vk0(i))*(rct1-two*rc*sm(i)+sm(i)**2)/rct2
187 ELSEIF (sm(i) >rok)THEN
188 vk(i)=vk0(i)
189 ELSEIF (sm(i) > rob(i))THEN ! Cape
190 vk(i)=vk0(i)*(one- ((sm(i)-rok)/(rob(i)-rok))**2)
191 ELSE
192 vk(i)=zero
193 ENDIF
194c CRITERE
195 IF (ibug == 0) THEN
196 aj3 = s1(i)*s2(i)*s3(i)
197 . - s1(i)*s5(i)*s5(i)-s2(i)*s6(i)*s6(i)-s3(i)*s4(i)*s4(i)
198 . + two*s4(i)*s5(i)*s6(i)
199 ELSE
200 aj3 = s1(i)*s2(i)*s3(i)
201 . - s1(i)*s5(i)*s5(i)-s2(i)*s6(i)*s6(i)-s3(i)*s4(i)*s4(i)
202 ENDIF
203 cs3t= half * aj3*(three/(half*max(r2,em20)))**three_half
204 cs3t= min(one,cs3t)
205 cs3t= max(-one,cs3t)
206 bb = half*((one - cs3t)*bc+(one +cs3t)*bt)
207 df = sqrt(bb*bb+ max(-aa*sm(i)+ac,em9))
208 rf = (-bb+df)/aa
209 rf2 = rf**2
210 aj2 = half*r2
211 ajj = sqrt(aj2)
212C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
213C NOW COMPUTES ALL TERMS TO BUILD UP PLASTIC MATRIX
214C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
215C D I M E N S I O N S RF:CONT TO:CONT PHI:ADIM B0:ADIM B1:1/CONT
216C B2:1/CONT**2 DFDTO:ADIM DRFDSM:ADIM DKDSM:1/CONT
217C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
218C PRINCIPALES DERIVEES
219 sm(i) = max(rob(i),sm(i))
220 IF(sm(i) >= rt-tol)THEN
221 dkdsm =zero
222 ELSEIF(sm(i)>rc)THEN
223 dkdsm = two*(one -vk0(i))*(sm(i)-rc)/rct2
224 ELSEIF(sm(i)>rok)THEN
225 dkdsm = zero
226 ELSE
227 dkdsm = -two*vk0(i)*(sm(i)-rok)/(rob(i)-rok)**2
228 ENDIF
229 drfdsm= -half / df
230 b0 = -third * vk(i) *drfdsm - third * rf * dkdsm
231 IF(ajj>em3*fc)THEN
232 drf3 = half* (-one + bb/df) * (bt-bc)/aa
233 b1 = half*sqr2/max(ajj,em20)
234 . + vk(i)*drf3*fourth*aj3*(three/max(aj2,em20))**twop5
235 b2 =-vk(i)*drf3*half*(three/max(aj2,em20))**three_half
236 ELSE
237 b1=zero
238 b2=zero
239 ENDIF
240C DF/DSIG
241 ts1 = s1(i)**2 + s4(i)**2 + s6(i)**2 - two*third*aj2
242 ts2 = s2(i)**2 + s4(i)**2 + s5(i)**2 - two*third*aj2
243 ts3 = s3(i)**2 + s5(i)**2 + s6(i)**2 - two*third*aj2
244 IF (ibug == 0) THEN
245 ts4 = two* (s5(i)*s6(i)-s4(i)*s3(i))
246 ts5 = two* (s6(i)*s4(i)-s5(i)*s1(i))
247 ts6 = two* (s4(i)*s5(i)-s6(i)*s2(i))
248 ELSE
249 ts4 = -two* s4(i) * s3(i)
250 ts5 = -two* s5(i) * s1(i)
251 ts6 = -two* s6(i) * s2(i)
252 ENDIF
253 dfs1=b0+b1*s1(i)+b2*ts1
254 dfs2=b0+b1*s2(i)+b2*ts2
255 dfs3=b0+b1*s3(i)+b2*ts3
256 dfs4=two*b1*s4(i)+b2*ts4
257 dfs5=two*b1*s5(i)+b2*ts5
258 dfs6=two*b1*s6(i)+b2*ts6
259C DG/DSIG
260C PLASTICITE VOLUMIQUE COMPACTANCE/DILATANCE
261 IF (vk(i) > vky) THEN
262 alpha = (one-vk(i))*ali0+(vk(i)-vky)*alf0
263 alpha = alpha/(one-vky)
264 ELSE
265 alpha = ali0
266 ENDIF
267 IF (icap == 1) THEN
268 IF (b0 < -em02) alpha = min(alpha,-zep4/b0)
269 IF (b0 > em02) alpha = max(alpha,-zep4/b0)
270 ENDIF
271 IF (eint(i)<=zero) alpha=zero
272c AJOUT DILATANCE PLASTIQUE MAX fp 09/16
273 IF (rho(i) < rho0) alpha = zero
274 IF (ajj > em3*fc) THEN
275 dgs1=alpha+s1(i)/(two*ajj)
276 dgs2=alpha+s2(i)/(two*ajj)
277 dgs3=alpha+s3(i)/(two*ajj)
278 dgs4=s4(i)/ajj
279 dgs5=s5(i)/ajj
280 dgs6=s6(i)/ajj
281 ELSE
282 IF (icap == 1) THEN
283 dgs1=alpha
284 dgs2=alpha
285 dgs3=alpha
286 ELSE
287 dgs1=-one
288 dgs2=-one
289 dgs3=-one
290 ENDIF
291 dgs4=zero
292 dgs5=zero
293 dgs6=zero
294 ENDIF
295c ECROUISSAGE
296 hpv = hv0*exp( min(fifty,(rob(i)-ro0)*expo))
297 IF(sm(i)>rok0)THEN
298 hp=hbp
299 ELSE
300 hp=hpv
301 ENDIF
302c
303 IF (icap == 1) THEN
304 phi = (alpha*three*sm(i)+ajj)
305 dfdto1=b0-sqrt(two_third)
306 dfdto2=three*b0
307 IF(dfdto1<=dfdto2)THEN
308 to=sqr3_2*vk(i)*rf
309 dfdto=dfdto1
310 ELSE
311 to=abs(sm(i))
312 dfdto=dfdto2
313 ENDIF
314 ecr=phi*hp*dfdto/to*(half-sign(half,vk(i)-one))
315 ELSE
316 dfdto1=b0-sqrt(two_third)
317 dfdto2=three*b0
318C IF(DFDTO2>=ZERO)THEN
319 IF(dfdto1<=dfdto2)THEN
320 to=sqr3_2*vk(i)*rf
321 phi = (alpha*three*sm(i)+ajj)/to
322 dfdto=dfdto1
323 ecr=phi*hp*dfdto*(half-sign(half,vk(i)-one))
324c write(IOUT,*)'BASE',SM(I),ROK,VK(I),DFDTO1,DFDTO2,ECR
325 ELSE
326 dfdro=-2*vk0(i)*rf*(sm(i)-rok)/(ro0-rok0)**2
327 ecr=dfdto2*dfdro*hpv
328c write(IOUT,*)'CAP',SM(I),ROK,VK(I),DFDTO1,DFDTO2,DFDRO,ECR
329 dgs1=dfs1
330 dgs2=dfs2
331 dgs3=dfs3
332 dgs4=dfs4
333 dgs5=dfs5
334 dgs6=dfs6
335 ENDIF
336 ENDIF
337c PRODUITS TENSORIELS
338 IF (icap == 1 .OR. (icap == 0 .AND. vk(i) > em05))THEN
339 h1a = cdam(i,1,1)*dfs1 + cdam(i,2,1)*dfs2 + cdam(i,3,1)*dfs3
340 h2a = cdam(i,1,2)*dfs1 + cdam(i,2,2)*dfs2 + cdam(i,3,2)*dfs3
341 h3a = cdam(i,1,3)*dfs1 + cdam(i,2,3)*dfs2 + cdam(i,3,3)*dfs3
342 h4a = c44(i)*dfs4
343 h5a = c55(i)*dfs5
344 h6a = c66(i)*dfs6
345C
346 h1n = cdam(i,1,1)*dgs1 + cdam(i,1,2)*dgs2 + cdam(i,1,3)*dgs3
347 h2n = cdam(i,2,1)*dgs1 + cdam(i,2,2)*dgs2 + cdam(i,2,3)*dgs3
348 h3n = cdam(i,3,1)*dgs1 + cdam(i,3,2)*dgs2 + cdam(i,3,3)*dgs3
349 h4n = c44(i)*dgs4
350 h5n = c55(i)*dgs5
351 h6n = c66(i)*dgs6
352c DENOMINATEUR
353 hh= dfs1*h1n+dfs2*h2n+dfs3*h3n+dfs4*h4n+dfs5*h5n+dfs6*h6n
354 . - min(zero,ecr)
355 hh = sign(max(abs(hh),em20),hh)
356c WRITE(IOUT,*)'HH',HH,ECR
357C
358 cp11=-h1n*h1a/hh*scal1(i)
359 cp12=-h1n*h2a/hh*scal1(i)*scal2(i)
360 cp13=-h1n*h3a/hh*scal1(i)*scal3(i)
361 cp14=-h1n*h4a/hh*scal1(i)*scal2(i)
362 cp15=-h1n*h5a/hh*scal1(i)*scal2(i)*scal3(i)
363 cp16=-h1n*h6a/hh*scal1(i)*scal3(i)
364C
365 cp21=-h2n*h1a/hh*scal1(i)*scal2(i)
366 cp22=-h2n*h2a/hh*scal2(i)
367 cp23=-h2n*h3a/hh*scal3(i)*scal2(i)
368 cp24=-h2n*h4a/hh*scal2(i)*scal1(i)
369 cp25=-h2n*h5a/hh*scal2(i)*scal3(i)
370 cp26=-h2n*h6a/hh*scal2(i)*scal1(i)*scal3(i)
371C
372 cp31=-h3n*h1a/hh*scal3(i)*scal1(i)
373 cp32=-h3n*h2a/hh*scal3(i)*scal2(i)
374 cp33=-h3n*h3a/hh*scal3(i)
375 cp34=-h3n*h4a/hh*scal3(i)*scal1(i)*scal2(i)
376 cp35=-h3n*h5a/hh*scal3(i)*scal2(i)
377 cp36=-h3n*h6a/hh*scal3(i)*scal1(i)
378C
379 cp41=-h4n*h1a/hh*scal1(i)
380 cp42=-h4n*h2a/hh*scal2(i)
381 cp43=-h4n*h3a/hh*scal3(i)
382 cp44=-h4n*h4a/hh*scal1(i)*scal2(i)
383 cp45=-h4n*h5a/hh*scal1(i)*scal2(i)*scal3(i)
384 cp46=-h4n*h6a/hh*scal1(i)*scal2(i)*scal3(i)
385C
386 cp51=-h5n*h1a/hh*scal1(i)
387 cp52=-h5n*h2a/hh*scal2(i)
388 cp53=-h5n*h3a/hh*scal3(i)
389 cp54=-h5n*h4a/hh*scal1(i)*scal2(i)*scal3(i)
390 cp55=-h5n*h5a/hh*scal2(i)*scal3(i)
391 cp56=-h5n*h6a/hh*scal1(i)*scal2(i)*scal3(i)
392C
393 cp61=-h6n*h1a/hh*scal1(i)
394 cp62=-h6n*h2a/hh*scal2(i)
395 cp63=-h6n*h3a/hh*scal3(i)
396 cp64=-h6n*h4a/hh*scal1(i)*scal2(i)*scal3(i)
397 cp65=-h6n*h5a/hh*scal1(i)*scal2(i)*scal3(i)
398 cp66=-h6n*h6a/hh*scal1(i)*scal3(i)
399C . . . . . . . . . . . . . . . . . . . . . . . . . . . .
400C NEW STRESS
401C DSIGI=SCLE2*(CIJ+CPIJ)EJ
402C . . . . . . . . . . . . . . . . . . . . . . . . . . . .
403 cp11=cdam(i,1,1)+cp11
404 cp12=cdam(i,1,2)+cp12
405 cp13=cdam(i,1,3)+cp13
406 cp21=cdam(i,2,1)+cp21
407 cp22=cdam(i,2,2)+cp22
408 cp23=cdam(i,2,3)+cp23
409 cp31=cdam(i,3,1)+cp31
410 cp32=cdam(i,3,2)+cp32
411 cp33=cdam(i,3,3)+cp33
412 cp44=c44(i)+cp44
413 cp55=c55(i)+cp55
414 cp66=c66(i)+cp66
415 sig(i,1)=sig(i,1)+cp11*e1(i)+cp12*e2(i)+cp13*e3(i)
416 . +cp14*e4(i)+cp15*e5(i)+cp16*e6(i)
417 sig(i,2)=sig(i,2)+cp21*e1(i)+cp22*e2(i)+cp23*e3(i)
418 . +cp24*e4(i)+cp25*e5(i)+cp26*e6(i)
419 sig(i,3)=sig(i,3)+cp31*e1(i)+cp32*e2(i)+cp33*e3(i)
420 . +cp34*e4(i)+cp35*e5(i)+cp36*e6(i)
421 sig(i,4)=sig(i,4)+cp41*e1(i)+cp42*e2(i)+cp43*e3(i)
422 . +cp44*e4(i)+cp45*e5(i)+cp46*e6(i)
423 sig(i,5)=sig(i,5)+cp51*e1(i)+cp52*e2(i)+cp53*e3(i)
424 . +cp54*e4(i)+cp55*e5(i)+cp56*e6(i)
425 sig(i,6)=sig(i,6)+cp61*e1(i)+cp62*e2(i)+cp63*e3(i)
426 . +cp64*e4(i)+cp65*e5(i)+cp66*e6(i)
427 ELSE
428c PURE TRIAXIAL TO SOLVE UNDETERMINATION
429 dp = bulk*hpv/(bulk+hpv)*(e1(i)+e2(i)+e3(i))
430 sig(i,1)=sig(i,1)+dp
431 sig(i,2)=sig(i,2)+dp
432 sig(i,3)=sig(i,3)+dp
433 ENDIF
434C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
435C RECALCUL DES DEFORMATION ELASTIQUES POUR REOUVERTURE DES FISSURES
436C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
437 c11 = one/de1(i)/young
438 c12 = -nu*scal1(i)*scal2(i)/young
439 c13 = -nu*scal1(i)*scal3(i)/young
440 crak(i,1) = c11*sig(i,1)+c12*sig(i,2)+c13*sig(i,3)
441 c22 = one/de2(i)/young
442 c23 = -nu*scal2(i)*scal3(i)/young
443 crak(i,2) = c12*sig(i,1)+c22*sig(i,2)+c23*sig(i,3)
444 c33 = one/de3(i)/young
445 crak(i,3) = c13*sig(i,1)+c23*sig(i,2)+c33*sig(i,3)
446C
447 de1(i)=one - max( zero , sign(dam(i,1),crak(i,1)) )
448 de2(i)=one - max( zero , sign(dam(i,2),crak(i,2)) )
449 de3(i)=one - max( zero , sign(dam(i,3),crak(i,3)) )
450 scal1(i)=half+sign(half,de1(i)-one)
451 scal2(i)=half+sign(half,de2(i)-one)
452 scal3(i)=half+sign(half,de3(i)-one)
453 de4(i)=scal1(i)*scal2(i)
454 de5(i)=scal2(i)*scal3(i)
455 de6(i)=scal3(i)*scal1(i)
456c--------------------------
457 den = one - nu**2 *(de4(i) + de5(i) + de6(i)
458 . + two*nu*scal1(i)*scal2(i)*scal3(i))
459C
460 cdam(i,1,1) = young*de1(i)*(one-nu**2*scal2(i)*scal3(i))/den
461 cdam(i,2,2) = young*de2(i)*(one -nu**2*scal1(i)*scal3(i))/den
462 cdam(i,3,3) = young*de3(i)*(one -nu**2*scal2(i)*scal1(i))/den
463 cdam(i,1,2) = nu*young*scal1(i)*scal2(i) *(one+nu*scal3(i))/den
464 cdam(i,2,3) = nu*young*scal2(i)*scal3(i) *(one+nu*scal1(i))/den
465 cdam(i,1,3) = nu*young*scal1(i)*scal3(i) *(one+nu*scal2(i))/den
466 cdam(i,2,1) = cdam(i,1,2)
467 cdam(i,3,1) = cdam(i,1,3)
468 cdam(i,3,2) = cdam(i,2,3)
469 sig(i,1) = cdam(i,1,1)*crak(i,1)+cdam(i,1,2)*crak(i,2)+cdam(i,1,3)*crak(i,3)
470 sig(i,2) = cdam(i,2,1)*crak(i,1)+cdam(i,2,2)*crak(i,2)+cdam(i,2,3)*crak(i,3)
471 sig(i,3) = cdam(i,3,1)*crak(i,1)+cdam(i,3,2)*crak(i,2)+cdam(i,3,3)*crak(i,3)
472 sig(i,4) = sig(i,4)*de4(i)
473 sig(i,5) = sig(i,5)*de5(i)
474 sig(i,6) = sig(i,6)*de6(i)
475c--------------------------
476 s1(i) = sig(i,1) * scal1(i)
477 s2(i) = sig(i,2) * scal2(i)
478 s3(i) = sig(i,3) * scal3(i)
479 s4(i) = sig(i,4) * de4(i)
480 s5(i) = sig(i,5) * de5(i)
481 s6(i) = sig(i,6) * de6(i)
482 sm(i) = third * (s1(i)+s2(i)+s3(i))
483C
484 s1(i)= s1(i)-sm(i)
485 s2(i)= s2(i)-sm(i)
486 s3(i)= s3(i)-sm(i)
487C
488C . . . . . . . . . . . . . . . . . . . . . . . . . . . .
489C CALCUL DES PARAMETRES D' ECROUISSAGE COURANTS VK0,ROB(I)
490C CORRECTIONS EVENTUELLES
491C . . . . . . . . . . . . . . . . . . . . . . . . . . . .
492 vkold=vk(i)
493 IF (sm(i)>ac/aa) sm(i) =
494 . sm(i) -three*(sm(i)-ac/aa)/(scal1(i)+scal2(i)+scal3(i))
495C
496 r2 = s1(i)**2+s2(i)**2+s3(i)**2
497 . + two*s4(i)**2+two*s5(i)**2+two*s6(i)**2
498 rr = sqrt(r2)
499 IF (ibug == 0) THEN
500 aj3 = s1(i)*s2(i)*s3(i)
501 . - s1(i)*s5(i)*s5(i)-s2(i)*s6(i)*s6(i)-s3(i)*s4(i)*s4(i)
502 . + two*s4(i)*s5(i)*s6(i)
503 ELSE
504 aj3 = s1(i)*s2(i)*s3(i)
505 . - s1(i)*s5(i)*s5(i)-s2(i)*s6(i)*s6(i)-s3(i)*s4(i)*s4(i)
506 ENDIF
507 cs3t= half * aj3*(three/(half*max(r2,em20)))**three_half
508 cs3t= min(one,cs3t)
509 cs3t= max(-one,cs3t)
510 bb = half*((one-cs3t)*bc+(one+cs3t)*bt)
511 df = sqrt(bb*bb+ max(-aa*sm(i)+ac,zero))
512 rf = (-bb+df)/aa
513 vk(i) = rr/max(rf,em20)
514 vkmax=one
515C
516c CORRECTION SI EN DEHORS DU CRITERE DE RUPTURE
517 IF (vk(i) > vkmax)THEN
518 fac = vkmax/vk(i)
519 IF(scal1(i) > zep9) sig(i,1) = s1(i)*fac+sm(i)
520 IF(scal2(i) > zep9) sig(i,2) = s2(i)*fac+sm(i)
521 IF(scal3(i) > zep9) sig(i,3) = s3(i)*fac+sm(i)
522 IF(scal1(i)*scal2(i) > zep9)sig(i,4) = s4(i)*fac
523 IF(scal2(i)*scal3(i) > zep9)sig(i,5) = s5(i)*fac
524 IF(scal3(i)*scal1(i) > zep9)sig(i,6) = s6(i)*fac
525 vk(i)=vkmax
526 c11 = one/de1(i)/young
527 c12 = -nu*scal1(i)*scal2(i)/young
528 c13 = -nu*scal1(i)*scal3(i)/young
529 crak(i,1) = c11*sig(i,1)+c12*sig(i,2)+c13*sig(i,3)
530 c22 = one/de2(i)/young
531 c23 = -nu*scal2(i)*scal3(i)/young
532 crak(i,2) = c12*sig(i,1)+c22*sig(i,2)+c23*sig(i,3)
533 c33 = one/de3(i)/young
534 crak(i,3) = c13*sig(i,1)+c23*sig(i,2)+c33*sig(i,3)
535C
536 de1(i)=one - max( zero , sign(dam(i,1),crak(i,1)) )
537 de2(i)=one - max( zero , sign(dam(i,2),crak(i,2)) )
538 de3(i)=one - max( zero , sign(dam(i,3),crak(i,3)) )
539 scal1(i)=half + sign(half,de1(i)-one)
540 scal2(i)=half + sign(half,de2(i)-one)
541 scal3(i)=half + sign(half,de3(i)-one)
542 de4(i)=scal1(i)*scal2(i)
543 de5(i)=scal2(i)*scal3(i)
544 de6(i)=scal3(i)*scal1(i)
545c-------------------------------
546 den = one - nu**2 *(de4(i) + de5(i) + de6(i)
547 . + two*nu*scal1(i)*scal2(i)*scal3(i))
548C
549 cdam(i,1,1) = young*de1(i)*(one -nu**2*scal2(i)*scal3(i))/den
550 cdam(i,2,2) = young*de2(i)*(one -nu**2*scal1(i)*scal3(i))/den
551 cdam(i,3,3) = young*de3(i)*(one -nu**2*scal2(i)*scal1(i))/den
552 cdam(i,1,2) = nu*young*scal1(i)*scal2(i) *(one +nu*scal3(i))/den
553 cdam(i,2,3) = nu*young*scal2(i)*scal3(i) *(one +nu*scal1(i))/den
554 cdam(i,1,3) = nu*young*scal1(i)*scal3(i) *(one +nu*scal2(i))/den
555 cdam(i,2,1) = cdam(i,1,2)
556 cdam(i,3,1) = cdam(i,1,3)
557 cdam(i,3,2) = cdam(i,2,3)
558 sig(i,1)=cdam(i,1,1)*crak(i,1)+cdam(i,1,2)*crak(i,2)+cdam(i,1,3)*crak(i,3)
559 sig(i,2)=cdam(i,2,1)*crak(i,1)+cdam(i,2,2)*crak(i,2)+cdam(i,2,3)*crak(i,3)
560 sig(i,3)=cdam(i,3,1)*crak(i,1)+cdam(i,3,2)*crak(i,2)+cdam(i,3,3)*crak(i,3)
561 sig(i,4)=de4(i)*sig(i,4)
562 sig(i,5)=de5(i)*sig(i,5)
563 sig(i,6)=de6(i)*sig(i,6)
564 ENDIF
565C
566 ro=rob(i)
567 IF (sm(i) >= rt-tol)THEN
568 vk(i)=one
569 ELSEIF(sm(i) > rc) THEN
570 div= min(-em20,rct1- two*rc*sm(i)+sm(i)*sm(i) )
571 vk(i)=one +(one - vk(i))*rct2/div
572 ELSEIF(sm(i) > rok) THEN
573 vk(i)=vk(i)
574 ELSE
575 dvk=vk(i)-vk0(i)*(one -(( max(sm(i),rob(i))-rok)/(ro0-rok0))**2)
576 vkk=vk0(i)+ max(dvk,zero)
577 vkk= min(vkk,one)
578 ro=sm(i)+(one -sqrt(one -vk(i)/vkk))*(ro0-rok0)
579 ENDIF
580 rob(i)= min(ro,rob(i))
581 vk(i) = min(vk(i),one)
582 vk0(i)= max(vk(i),vk0(i))
583c-----------
584c FIN DES ITERATIONS
585 ENDDO ! ITER
586c-----------
587 ENDDO ! NINDX
588c-----------
589 RETURN
590 END
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
plas24
subroutine plas24(nel, nindx, indx, ngl, pm, sig, dam, crak, rho, eint, vk0, vk, rob, cdam, e1, e2, e3, e4, e5, e6, s1, s2, s3, s4, s5, s6, scal1, scal2, scal3, scle2)
Definition plas24.F:34