37
38
39
40#include "implicit_f.inc"
41
42
43
44#include "mvsiz_p.inc"
45
46
47
48#include "com01_c.inc"
49
50
51
52
53 INTEGER NEL,NUPARAM,NUVAR,INLOC
55 . time,timestep(nel),uparam(nuparam),
56 . rho0(nel),thkly(nel),pla(nel),
57 . depsxx(nel),depsyy(nel),
58 . depsxy(nel),depsyz(nel),depszx(nel),
59 . sigoxx(nel),sigoyy(nel),
60 . sigoxy(nel),sigoyz(nel),sigozx(nel),
61 . gs(*),hardm(nel),seq(nel),dmg(nel),
62 . dplanl(nel)
63 my_real,
DIMENSION(NEL),
INTENT(IN) :: loff
64
65
66
68 . signxx(nel),signyy(nel),
69 . signxy(nel),signyz(nel),signzx(nel),
70 . soundsp(nel),etse(nel),dpla(nel)
71
72
73
75 . uvar(nel,nuvar), off(nel),thk(nel),yld(nel)
76
77
78
79 INTEGER I,II,J,NMAX,NINDX,INDEX(MVSIZ),NITER,ITER
81 . e,nu,g,a11,a12,
82 . sigy,eps0,nexp,
83 . ff,gg,hh,nn,
84 . c1,c2,c3,mexp,dc
86 . eta,cos3theta,theta,f1,f2,f3,cc(mvsiz),
87 . beta(mvsiz),dam(mvsiz),dezz(mvsiz),
88 . epsf,p,sighl(mvsiz),h(mvsiz),svm(mvsiz),
89 . normxx,normyy,normxy,dpxx(mvsiz),dpyy(mvsiz),
90 . dpzz(mvsiz),dpxy(mvsiz),deelzz(mvsiz),
91 . dpla_dlam(mvsiz),dlam,sig_dfdsig,dfdsig2,
92 . phi(mvsiz),dphi_dlam(mvsiz),ddep
93
94
95
96
97
98 e = uparam(1)
99 nu = uparam(2) ! poisson's ratio
100 G = UPARAM(3) ! Shear modulus
101 A11 = UPARAM(5) ! Diag. component of plane stress elastic matrix
102 A12 = UPARAM(6) ! Non diag component of plane stress elastic matrix
103 ! Hardening parameters
104 SIGY = UPARAM(11) ! Initial yield stress
105 EPS0 = UPARAM(12) ! Initial plastic strain (> 0)
106 NEXP = UPARAM(13) ! Hardening exponent
107 ! Hill yield criterion parameters
108 FF = UPARAM(14) ! First Hill coefficient
109 GG = UPARAM(15) ! Second Hill coefficient
110 HH = UPARAM(16) ! Third Hill coefficient
111 NN = UPARAM(17) ! Fourth Hill coefficient
112 ! Modified Mohr-Coulomb yield criterion parameters
113 C1 = UPARAM(20) ! First failure parameter
114 C2 = UPARAM(21) ! Second failure parameter
115 C3 = UPARAM(22) ! Third failure parameter
116 MEXP = UPARAM(23) ! Damage exponent
117 DC = UPARAM(24) ! Critical damage value
118
119 ! Initial value
120 IF (ISIGI == 0) THEN
121 IF (TIME == ZERO) THEN
122 DO I=1,NEL
123 UVAR(I,1) = ZERO
124 ENDDO
125 ENDIF
126 ENDIF
127
128 ! Recovering internal variables
129 DO I=1,NEL
130 ! Checking deletion flag value
131 IF (OFF(I) < ONE) OFF(I) = FOUR_OVER_5*OFF(I)
132 IF (OFF(I) < EM01) OFF(I) = ZERO
133 ! Hourglass coefficient
134 ETSE(I) = ONE
135 H(I) = ZERO
136 ! Plastic strain increment
137 DPLA(I) = ZERO
138 DPXX(I) = ZERO
139 DPYY(I) = ZERO
140 DPZZ(I) = ZERO
141 DPXY(I) = ZERO
142 ! Damage variable
143 DAM(I) = UVAR(I,1)
144 ENDDO
145
146 ! Computing yield stress and checking damage criteria
147 DO I = 1,NEL
148 ! Computation of the BETA factor
149 BETA(I) = ONE
150 IF (OFF(I) == ONE) THEN
151.AND. IF ((DAM(I) <= DC) (DAM(I) >= ONE)) THEN
152 BETA(I) = ONE/(MAX(DC - ONE,EM20))
153 BETA(I) = (DC - DAM(I))*BETA(I)
154 BETA(I) = BETA(I)**MEXP
155 ENDIF
156 ENDIF
157 ! Computation of the yield stress
158 CC(I) = PLA(I) + EPS0
159 YLD(I) = BETA(I)*SIGY*EXP(NEXP*LOG(CC(I)))
160 YLD(I) = MAX(YLD(I), EM10)
161 ENDDO
162
163 !========================================================================
164 ! - COMPUTATION OF TRIAL VALUES
165 !========================================================================
166 DO I = 1,NEL
167
168 ! Computation of the trial stress tensor
169 SIGNXX(I) = SIGOXX(I) + A11*DEPSXX(I)+A12*DEPSYY(I)
170 SIGNYY(I) = SIGOYY(I) + A12*DEPSXX(I)+A11*DEPSYY(I)
171 SIGNXY(I) = SIGOXY(I) + G*DEPSXY(I)
172 SIGNYZ(I) = SIGOYZ(I) + GS(I)*DEPSYZ(I)
173 SIGNZX(I) = SIGOZX(I) + GS(I)*DEPSZX(I)
174
175 ! Hill equivalent stress
176 SIGHL(I) = (FF+HH)*SIGNYY(I)**2 + (GG+HH)*SIGNXX(I)**2
177 . - TWO*HH*SIGNXX(I)*SIGNYY(I) + TWO*NN*SIGNXY(I)**2
178 SIGHL(I) = SQRT(MAX(ZERO,SIGHL(I)))
179
180 ENDDO
181
182 !========================================================================
183 ! - COMPUTATION OF YIELD FONCTION
184 !========================================================================
185 PHI(1:NEL) = SIGHL(1:NEL) - YLD(1:NEL)
186 ! Checking plastic behavior for all elements
187 NINDX = 0
188 INDEX = 0
189 DO I=1,NEL
190.AND. IF ((PHI(I)>ZERO)(OFF(I) == ONE)) THEN
191 NINDX = NINDX+1
192 INDEX(NINDX) = I
193 ENDIF
194 ENDDO
195
196 !====================================================================
197 ! - PLASTIC CORRECTION WITH CUTTING PLANE ALGORITHM (NEWTON ITERATION)
198 !====================================================================
199
200 ! Number of maximum Newton iterations
201 NITER = 3
202
203 ! Loop over the iterations
204 DO ITER = 1, NITER
205#include "vectorize.inc"
206 ! Loop over yielding elements
207 DO II=1,NINDX
208 I = INDEX(II)
209
210 ! Note : in this part, the purpose is to compute for each iteration
211 ! a plastic multiplier allowing to update internal variables to satisfy
212 ! the consistency condition using the backward Euler implicit method
213 ! with a Newton-Raphson iterative procedure
214 ! Its expression at each iteration is : DLAMBDA = - PHI/DPHI_DLAMBDA
215 ! -> PHI : current value of yield function (known)
216 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
217 ! into account of internal variables kinetic :
218 ! plasticity, temperature and damage (to compute)
219
220 ! 1 - Computation of DPHI_DSIG the normal to the yield surface
221 !-------------------------------------------------------------
222 NORMXX = (GG*SIGNXX(I) + HH*(SIGNXX(I)-SIGNYY(I)))/(MAX(SIGHL(I),EM20))
223 NORMYY = (FF*SIGNYY(I) + HH*(SIGNYY(I)-SIGNXX(I)))/(MAX(SIGHL(I),EM20))
224 NORMXY = TWO*NN*SIGNXY(I)/(MAX(SIGHL(I),EM20))
225
226 ! 2 - Computation of DPHI_DLAMBDA
227 !---------------------------------------------------------
228
229 ! a) Derivative with respect stress increments tensor DSIG
230 ! --------------------------------------------------------
231 DFDSIG2 = NORMXX * (A11*NORMXX + A12*NORMYY)
232 . + NORMYY * (A11*NORMYY + A12*NORMXX)
233 . + NORMXY * NORMXY * G
234
235 ! b) Derivatives with respect to plastic strain P
236 ! ------------------------------------------------
237
238 ! i) Derivative of the yield stress with respect to plastic strain dYLD / dPLA
239 ! ----------------------------------------------------------------------------
240 H(I) = BETA(I)*NEXP*SIGY*EXP((NEXP-1)*LOG(CC(I)))
241
242 ! ii) Derivative of dPLA with respect to DLAM
243 ! -------------------------------------------
244 SIG_DFDSIG = SIGNXX(I) * NORMXX
245 . + SIGNYY(I) * NORMYY
246 . + SIGNXY(I) * NORMXY
247 DPLA_DLAM(I) = SIG_DFDSIG / MAX(YLD(I),EM20)
248
249 ! 3 - Computation of plastic multiplier and variables update
250 !----------------------------------------------------------
251
252 ! Derivative of PHI with respect to DLAM
253 DPHI_DLAM(I) = - DFDSIG2 - H(I)*DPLA_DLAM(I)
254 DPHI_DLAM(I) = SIGN(MAX(ABS(DPHI_DLAM(I)),EM20) ,DPHI_DLAM(I))
255
256 ! Computation of the plastic multiplier
257 DLAM = -PHI(I)/DPHI_DLAM(I)
258
259 ! Plastic strains tensor update
260 DPXX(I) = DLAM * NORMXX
261 DPYY(I) = DLAM * NORMYY
262 DPXY(I) = DLAM * NORMXY
263
264 ! Elasto-plastic stresses update
265 SIGNXX(I) = SIGNXX(I) - (A11*DPXX(I) + A12*DPYY(I))
266 SIGNYY(I) = SIGNYY(I) - (A11*DPYY(I) + A12*DPXX(I))
267 SIGNXY(I) = SIGNXY(I) - DPXY(I)*G
268
269 ! Cumulated plastic strain and strain rate update
270 DDEP = DLAM*DPLA_DLAM(I)
271 DPLA(I) = MAX(ZERO, DPLA(I) + DDEP)
272 PLA(I) = MAX(PLA(I) + DDEP,ZERO)
273
274 ! Update the hardening yield stress
275 CC(I) = PLA(I) + EPS0
276 YLD(I) = BETA(I)*SIGY*EXP(NEXP*LOG(CC(I)))
277 YLD(I) = MAX(YLD(I),EM10)
278
279 ! Update Hill equivalent stress
280 SIGHL(I) = (FF+HH)*SIGNYY(I)**2 + (GG+HH)*SIGNXX(I)**2
281 . - TWO*HH*SIGNXX(I)*SIGNYY(I) + TWO*NN*SIGNXY(I)**2
282 SIGHL(I) = SQRT(MAX(ZERO,SIGHL(I)))
283
284 ! Update yield function value
285 PHI(I) = SIGHL(I) - YLD(I)
286
287 ! Transverse strain update
288 DPZZ(I) = DPZZ(I) - (DPXX(I)+DPYY(I))
289
290 ENDDO
291 ! End of the loop over the yielding elements
292 ENDDO
293 ! End of the loop over the iterations
294 !===================================================================
295 ! - END OF PLASTIC CORRECTION WITH CUTTING PLANE ALGORITHM
296 !===================================================================
297
298 ! Update and store new internal variables
299 DO I=1,NEL
300 ! Modified Mohr Criteria Failure Model
301 IF (OFF(I) == ONE) THEN
302 ! Pressure
303 P = THIRD*(SIGNXX(I) + SIGNYY(I))
304 ! Von Mises equivalent stress
305 SVM(I) = SIGNXX(I)**2 + SIGNYY(I)**2 - SIGNXX(I)*SIGNYY(I) +
306 . THREE*SIGNXY(I)**2
307 SVM(I) = SQRT(SVM(I))
308
309 ! Computation of triaxiality
310 ETA = P/MAX(SVM(I),EM20)
311 IF (ETA < -TWO_THIRD) ETA = -TWO_THIRD
312 IF (ETA > TWO_THIRD) ETA = TWO_THIRD
313 ! Computation of Lode angle
314 COS3THETA = -HALF*TWENTY7*ETA*(ETA**2 - THIRD)
315 IF (COS3THETA < -ONE ) COS3THETA = -ONE
316 IF (COS3THETA > ONE ) COS3THETA = ONE
317 THETA = ONE - TWO*ACOS(COS3THETA)/PI
318
319 ! Computation of the failure factor
320 F1 = COS(THETA*PI/SIX)
321 F2 = SIN(THETA*PI/SIX)
322 F3 = C3 + (SQRT(THREE)/(TWO - SQRT(THREE)))*(ONE - C3)*(ONE/MAX(F1,EM20) - ONE)
323
324 ! Computation of the failure plastic strain
325 EPSF = (SIGY/MAX(C2,EM20))*F3*(F1*SQRT(THIRD*(ONE + C1**2)) + C1*(ETA + F2*THIRD))
326 EPSF = MAX(EPSF,EM20)**(-ONE/NEXP)
327
328 ! Computation of the damage variable
329 IF (INLOC > 0) THEN
330 DAM(I) = DAM(I) + MAX(DPLANL(I),ZERO)/MAX(EM20,EPSF)
331 ELSE
332 DAM(I) = DAM(I) + DPLA(I)/MAX(EM20,EPSF)
333 ENDIF
334 IF (DAM(I) >= DC) THEN
335 DAM(I) = DC
336 OFF(I) = FOUR_OVER_5
337 ENDIF
338
339 ! Store the new value of damage
340 UVAR(I,1) = DAM(I)
341 ENDIF
342 ! Equivalent stress
343 SEQ(I) = SIGHL(I)
344 ! Normalized damage
345 DMG(I) = DAM(I)/DC
346 ! Sound-speed
347 SOUNDSP(I) = SQRT(A11/RHO0(I))
348 ! Coefficient for hourglass
349 IF (DPLA(I) > ZERO) THEN
350 ETSE(I) = H(I) / (H(I) + E)
351 HARDM(I) = H(I)
352 ELSE
353 ETSE(I) = ONE
354 HARDM(I) = ZERO
355 ENDIF
356 ! Computation of the thickness variation
357 DEELZZ(I) = -NU*(SIGNXX(I)-SIGOXX(I)+SIGNYY(I)-SIGOYY(I))/E
358 IF (INLOC > 0) THEN
359 IF (LOFF(I) == ONE) THEN
360 NORMXX = (GG*SIGNXX(I) + HH*(SIGNXX(I)-SIGNYY(I)))/(MAX(SIGHL(I),EM20))
361 NORMYY = (FF*SIGNYY(I) + HH*(SIGNYY(I)-SIGNXX(I)))/(MAX(SIGHL(I),EM20))
362 DEZZ(I) = DEELZZ(I) - MAX(DPLANL(I),ZERO)*(NORMXX+NORMYY)
363 ENDIF
364 ELSE
365 DEZZ(I) = DEELZZ(I) + DPZZ(I)
366 ENDIF
367 THK(I) = THK(I) + DEZZ(I)*THKLY(I)*OFF(I)
368 ENDDO
369
1.15.0