OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat107_newton.F File Reference
#include "implicit_f.inc"
#include "com04_c.inc"
#include "vectorize.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine mat107_newton (nel, ngl, nuparam, nuvar, time, timestep, uparam, uvar, jthe, off, rho0, rho, pla, dpla, epsd, soundsp, depsxx, depsyy, depszz, depsxy, depsyz, depszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, sigy, et, nvartmp, numtabl, vartmp, itable, table)

Function/Subroutine Documentation

◆ mat107_newton()

subroutine mat107_newton ( integer nel,
integer, dimension(nel), intent(in) ngl,
integer nuparam,
integer nuvar,
time,
timestep,
intent(in) uparam,
intent(inout) uvar,
integer jthe,
intent(inout) off,
intent(in) rho0,
intent(in) rho,
intent(inout) pla,
intent(inout) dpla,
intent(inout) epsd,
intent(out) soundsp,
intent(in) depsxx,
intent(in) depsyy,
intent(in) depszz,
intent(in) depsxy,
intent(in) depsyz,
intent(in) depszx,
intent(in) sigoxx,
intent(in) sigoyy,
intent(in) sigozz,
intent(in) sigoxy,
intent(in) sigoyz,
intent(in) sigozx,
intent(out) signxx,
intent(out) signyy,
intent(out) signzz,
intent(out) signxy,
intent(out) signyz,
intent(out) signzx,
intent(out) sigy,
intent(out) et,
integer nvartmp,
integer numtabl,
integer, dimension(nel,nvartmp), intent(inout) vartmp,
integer, dimension(numtabl) itable,
type(ttable), dimension(ntable) table )

Definition at line 33 of file mat107_newton.F.

42 !=======================================================================
43 ! Modules
44 !=======================================================================
45 USE table_mod
46 USE interface_table_mod
47 !=======================================================================
48 ! Implicit types
49 !=======================================================================
50#include "implicit_f.inc"
51 !=======================================================================
52 ! Common
53 !=======================================================================
54#include "com04_c.inc"
55 !=======================================================================
56 ! Dummy arguments
57 !=======================================================================
58 INTEGER NEL,NUPARAM,NUVAR,JTHE,NUMTABL,ITABLE(NUMTABL),NVARTMP
59 INTEGER ,DIMENSION(NEL), INTENT(IN) :: NGL
60 my_real
61 . time,timestep
62 INTEGER, INTENT(INOUT) :: VARTMP(NEL,NVARTMP)
63 my_real,DIMENSION(NUPARAM), INTENT(IN) ::
64 . uparam
65 my_real,DIMENSION(NEL), INTENT(IN) ::
66 . rho,rho0,
67 . depsxx,depsyy,depszz,depsxy,depsyz,depszx,
68 . sigoxx,sigoyy,sigozz,sigoxy,sigoyz,sigozx
69c
70 my_real ,DIMENSION(NEL), INTENT(OUT) ::
71 . soundsp,sigy,et,
72 . signxx,signyy,signzz,signxy,signyz,signzx
73c
74 my_real ,DIMENSION(NEL), INTENT(INOUT) ::
75 . dpla,off
76 my_real ,DIMENSION(NEL,6),INTENT(INOUT) ::
77 . pla,epsd
78 my_real ,DIMENSION(NEL,NUVAR), INTENT(INOUT) ::
79 . uvar
80c
81 TYPE(TTABLE), DIMENSION(NTABLE) :: TABLE
82 !=======================================================================
83 ! Local Variables
84 !=======================================================================
85 INTEGER I,K,II,NINDX,INDEX(NEL),NITER,ITER,ITAB,ISMOOTH,IPOS(NEL,2)
86c
87 my_real
88 . young1,young2,young3,nu12,nu21,g12,a11,a12,a21,a22,c1,
89 . xi1,xi2,k1,k2,k3,k4,k5,k6,g1c,d1,d2,sigy1,
90 . cini1,s1,sigy2,cini2,s2,sigy1c,cini1c,s1c,
91 . sigy2c,cini2c,s2c,sigyt,cinit,st,g23,g31
92 my_real
93 . normsig,
94 . h,dpdt,dlam,ddep,depxx,depyy
95 my_real
96 . normxx,normyy,normxy,normzz,normyz,normzx,
97 . normpxx,normpyy,normpxy,normpzz,normpyz,normpzx,
98 . dphi_dlam,dphi,dfdsig2,dphi_dpla,dpxx,dpyy,dpxy,
99 . dphi_dr1,dphi_dr2,dphi_dr1c,dphi_dr2c,dphi_drt
100 my_real
101 . xscale1,yscale1,xscale1c,yscale1c,xscale2,yscale2,xscale2c,yscale2c,
102 . xscalet,yscalet,xvec(nel,2),asrate
103c
104 my_real, DIMENSION(NEL) ::
105 . alpha1,alpha2,alpha3,alpha4,alpha5,r1,r2,r1c,r2c,rt,
106 . dsigxx,dsigyy,dsigxy,dsigyz,dsigzx,hardp,phi,eta1,eta2,
107 . dr1_dp,dr2_dp,dr1c_dp,dr2c_dp,drt_dp,hardr,dpla2,dpla3,dpla4,
108 . dpla5,dpla6,beta1,beta2
109c
110 !=======================================================================
111 ! - INITIALISATION OF COMPUTATION ON TIME STEP
112 !=======================================================================
113 ! Recovering model parameters
114 ! Elastic parameters
115 young1 = uparam(1) ! Young modulus in direction 1 (MD)
116 young2 = uparam(2) ! Young modulus in direction 2 (CD)
117 young3 = uparam(3) ! Young modulus in direction 3 (Out-of-plane)
118 nu12 = uparam(4) ! Poisson's ratio in 12
119 nu21 = uparam(5) ! Poisson's ratio in 21
120 a11 = uparam(6) ! Component 11 of orthotropic shell elasticity matrix
121 a12 = uparam(7) ! Component 12 of orthotropic shell elasticity matrix
122 a21 = uparam(8) ! Component 21 of orthotropic shell elasticity matrix
123 a22 = uparam(9) ! Component 22 of orthotropic shell elasticity matrix
124 g12 = uparam(10) ! Shear modulus in 12
125 g23 = uparam(11) ! shear modulus in 23
126 g31 = uparam(12) ! Shear modulus in 31
127 itab = int(uparam(14)) ! Tabulated yield stress flag
128 xi1 = uparam(15) ! 1st coupling parameter
129 xi2 = uparam(16) ! 2nd coupling parameter
130 k1 = uparam(17) ! 1st plastic potential parameter
131 k2 = uparam(18) ! 2nd plastic potential parameter
132 k3 = uparam(19) ! 3rd plastic potential parameter
133 k4 = uparam(20) ! 4th plastic potential parameter
134 k5 = uparam(21) ! 5th plastic potential parameter
135 k6 = uparam(22) ! 6th plastic potential parameter
136 g1c = uparam(23) ! Correction factor for R1C
137 IF (itab == 0) THEN
138 d1 = uparam(24) ! 1st Shear yield stress parameter
139 d2 = uparam(25) ! 2nd Shear yield stress parameter
140 sigy1 = uparam(26) ! Initial yield stress in tension in direction 1 (MD)
141 cini1 = uparam(27) ! Yield stress intersection with ordinate axis in tension in direction 1 (MD)
142 s1 = uparam(28) ! Yield stress slope in tension in direction 1 (MD)
143 sigy2 = uparam(29) ! Initial yield stress in tension in direction 2 (CD)
144 cini2 = uparam(30) ! Yield stress intersection with ordinate axis in tension in direction 2 (CD)
145 s2 = uparam(31) ! Yield stress slope in tension in direction 2 (CD)
146 sigy1c = uparam(32) ! Initial yield stress in compression in direction 1 (MD)
147 cini1c = uparam(33) ! Yield stress intersection with ordinate axis in compression in direction 1 (MD)
148 s1c = uparam(34) ! Yield stress slope in compression in direction 1 (MD)
149 sigy2c = uparam(35) ! Initial yield stress in compression in direction 2 (CD)
150 cini2c = uparam(36) ! Yield stress intersection with ordinate axis in compression in direction 2 (CD)
151 s2c = uparam(37) ! yield stress slope in compression in direction 2 (cd)
152 sigyt = uparam(38) ! Initial yield stress in shear
153 cinit = uparam(39) ! Yield stress intersection with ordinate axis in shear
154 st = uparam(40) ! Yield stress slope in shear
155 ELSE
156 xscale1 = uparam(24)
157 yscale1 = uparam(25)
158 xscale2 = uparam(26)
159 yscale2 = uparam(27)
160 xscale1c= uparam(28)
161 yscale1c= uparam(29)
162 xscale2c= uparam(30)
163 yscale2c= uparam(31)
164 xscalet = uparam(32)
165 yscalet = uparam(33)
166 asrate = uparam(34)
167 asrate = (two*pi*asrate*timestep)/(two*pi*asrate*timestep + one)
168 ismooth = int(uparam(35))
169 ENDIF
170c
171 ! Recovering internal variables
172 DO i=1,nel
173 ! OFF parameter for element deletion
174 IF (off(i) < 0.1) off(i) = zero
175 IF (off(i) < one) off(i) = off(i)*four_over_5
176 ! Standard inputs
177 dpla(i) = zero ! Initialization of the "global" plastic strain increment
178 dpla2(i) = zero ! Initialization of the in-plane plastic strain increment
179 dpla3(i) = zero ! Initialization of the transverse shear plastic strain increment
180 dpla4(i) = zero ! Initialization of the in-plane plastic strain increment
181 dpla5(i) = zero ! Initialization of the transverse shear plastic strain increment
182 dpla6(i) = zero ! Initialization of the in-plane plastic strain increment
183 et(i) = one ! initialization of hourglass coefficient
184 hardp(i) = zero ! Initialization of hourglass coefficient
185 ENDDO
186c
187 ! Computation of the initial yield stress
188 IF (itab > 0) THEN
189 ! -> Tensile yield stress in direction 1 (MD)
190 xvec(1:nel,1) = pla(1:nel,2)
191 xvec(1:nel,2) = epsd(1:nel,2) * xscale1
192 ipos(1:nel,1) = vartmp(1:nel,1)
193 ipos(1:nel,2) = 1
194 CALL table2d_vinterp_log(table(itable(1)),ismooth,nel,nel,ipos,xvec,r1,dr1_dp,hardr)
195 r1(1:nel) = r1(1:nel) * yscale1
196 dr1_dp(1:nel) = dr1_dp(1:nel) * yscale1
197 vartmp(1:nel,1) = ipos(1:nel,1)
198 ! -> compressive yield stress in direction 1 (md)
199 xvec(1:nel,1) = pla(1:nel,3)
200 xvec(1:nel,2) = epsd(1:nel,3) * xscale2
201 ipos(1:nel,1) = vartmp(1:nel,2)
202 ipos(1:nel,2) = 1
203 CALL table2d_vinterp_log(table(itable(2)),ismooth,nel,nel,ipos,xvec,r2,dr2_dp,hardr)
204 r2(1:nel) = r2(1:nel) * yscale2
205 dr2_dp(1:nel) = dr2_dp(1:nel) * yscale2
206 vartmp(1:nel,2) = ipos(1:nel,1)
207 ! -> Positive shear yield stress
208 xvec(1:nel,1) = pla(1:nel,4)
209 xvec(1:nel,2) = epsd(1:nel,4) * xscale1c
210 ipos(1:nel,1) = vartmp(1:nel,3)
211 ipos(1:nel,2) = 1
212 CALL table2d_vinterp_log(table(itable(3)),ismooth,nel,nel,ipos,xvec,r1c,dr1c_dp,hardr)
213 r1c(1:nel) = r1c(1:nel) * yscale1c
214 dr1c_dp(1:nel)= dr1c_dp(1:nel) * yscale1c
215 vartmp(1:nel,3) = ipos(1:nel,1)
216 ! -> Tensile yield stress in direction 2 (CD)
217 xvec(1:nel,1) = pla(1:nel,5)
218 xvec(1:nel,2) = epsd(1:nel,5) * xscale2c
219 ipos(1:nel,1) = vartmp(1:nel,4)
220 ipos(1:nel,2) = 1
221 CALL table2d_vinterp_log(table(itable(4)),ismooth,nel,nel,ipos,xvec,r2c,dr2c_dp,hardr)
222 r2c(1:nel) = r2c(1:nel) * yscale2c
223 dr2c_dp(1:nel)= dr2c_dp(1:nel) * yscale2c
224 vartmp(1:nel,4) = ipos(1:nel,1)
225 ! Compressive yield stress in direction 2 (CD)
226 xvec(1:nel,1) = pla(1:nel,6)
227 xvec(1:nel,2) = epsd(1:nel,6) * xscalet
228 ipos(1:nel,1) = vartmp(1:nel,5)
229 ipos(1:nel,2) = 1
230 CALL table2d_vinterp_log(table(itable(5)),ismooth,nel,nel,ipos,xvec,rt,drt_dp,hardr)
231 rt(1:nel) = rt(1:nel) * yscalet
232 drt_dp(1:nel) = drt_dp(1:nel) * yscalet
233 vartmp(1:nel,5) = ipos(1:nel,1)
234 ENDIF
235c
236 !========================================================================
237 ! - COMPUTATION OF TRIAL VALUES
238 !========================================================================
239 DO i=1,nel
240c
241 ! Computation of the trial stress tensor
242 signxx(i) = sigoxx(i) + a11*depsxx(i) + a12*depsyy(i)
243 signyy(i) = sigoyy(i) + a21*depsxx(i) + a22*depsyy(i)
244 signzz(i) = sigozz(i) + young3*depszz(i)
245 signxy(i) = sigoxy(i) + depsxy(i)*g12
246 signyz(i) = sigoyz(i) + depsyz(i)*g23
247 signzx(i) = sigozx(i) + depszx(i)*g31
248c
249 ! Computation of trial alpha coefficients
250 alpha1(i) = zero
251 alpha2(i) = zero
252 alpha3(i) = zero
253 alpha4(i) = zero
254 alpha5(i) = zero
255 IF ((signxx(i) - xi1*signyy(i)) > zero) alpha1(i) = one
256 IF ((signyy(i) - xi2*signxx(i)) > zero) alpha2(i) = one
257 IF (-signxx(i) > zero) alpha3(i) = one
258 IF (-signyy(i) > zero) alpha4(i) = one
259 IF (abs(signxy(i)) > zero) alpha5(i) = one
260c
261 IF (itab == 0) THEN
262 ! Computation of yield stresses
263 ! Tensile yield stresses
264 ! - in direction 1 (MD)
265 r1(i) = sigy1 + (one/(cini1 + s1*pla(i,2)))*pla(i,2)
266 ! - in direction 2 (CD)
267 r2(i) = sigy2 + (one/(cini2 + s2*pla(i,3)))*pla(i,3)
268 ! Compressive yield stresses
269 ! - in direction 1 (MD) (including correction)
270 r1c(i) = sigy1c + (one/(cini1c + s1c*pla(i,4)))*pla(i,4)
271 r1c(i) = r1c(i)/sqrt(one - g1c)
272 ! - in direction 2 (CD)
273 r2c(i) = sigy2c + (one/(cini2c + s2c*pla(i,5)))*pla(i,5)
274 ! Shear yield stress
275 eta1(i) = one + alpha3(i)*alpha4(i)*d1
276 eta2(i) = one + alpha3(i)*alpha4(i)*d2
277 rt(i) = eta1(i)*sigyt + eta2(i)*(one/(cinit + st*pla(i,6)))*pla(i,6)
278 ENDIF
279c
280 ! Computation of BETA switching coefficient
281 normsig = sqrt(signxx(i)*signxx(i) + signyy(i)*signyy(i) + two*signxy(i)*signxy(i))
282 normsig = max(normsig,em20)
283 beta1(i) = zero
284 beta2(i) = zero
285 IF ((signxx(i)/normsig > em01).OR.(signyy(i)/normsig > em01)) beta1(i) = one
286 IF ((signxx(i)/normsig < -em01).OR.(signyy(i)/normsig < -em01)) beta2(i) = one
287 IF (((abs(signxx(i))/normsig)<em01).AND.((abs(signyy(i))/normsig)<em01)) THEN
288 beta1(i) = one
289 beta2(i) = one
290 ENDIF
291c
292 ! Computation of the yield function
293 phi(i) = alpha1(i)*((signxx(i) - xi1*signyy(i))/r1(i))**2 +
294 . alpha2(i)*((signyy(i) - xi2*signxx(i))/r2(i))**2 +
295 . alpha3(i)*((- signxx(i))/r1c(i))**2 +
296 . alpha4(i)*((- signyy(i))/r2c(i))**2 +
297 . alpha5(i)*(abs(signxy(i))/rt(i))**2 - one
298 ENDDO
299c
300 !========================================================================
301 ! - CHECKING THE PLASTIC BEHAVIOR
302 !========================================================================
303c
304 ! Checking plastic behavior for all elements
305 nindx = 0
306 DO i=1,nel
307 IF (phi(i) > zero .AND. off(i) == one) THEN
308 nindx=nindx+1
309 index(nindx)=i
310 ENDIF
311 ENDDO
312c
313 !============================================================
314 ! - PLASTIC CORRECTION WITH CUTTING PLANE (NEWTON RESOLUTION)
315 !============================================================
316c
317 ! Number of Newton iterations
318 niter = 3
319c
320 ! Loop over the iterations
321 DO iter = 1, niter
322#include "vectorize.inc"
323 ! Loop over yielding elements
324 DO ii=1,nindx
325 i = index(ii)
326c
327 ! Note : in this part, the purpose is to compute for each iteration
328 ! a plastic multiplier allowing to update internal variables to satisfy
329 ! the consistency condition using the cutting plane method
330 ! Its expression at each iteration is : DLAMBDA = - PHI/DPHI_DLAMBDA
331 ! -> PHI : current value of yield function (known)
332 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
333 ! into account of internal variables kinetic :
334 ! hardening parameters
335c
336 ! 1 - Computation of DPHI_DSIG the normal to the yield criterion
337 !-------------------------------------------------------------
338c
339 normxx = two*(alpha1(i)/(r1(i)**2))*(signxx(i)-xi1*signyy(i))
340 . - two*(alpha2(i)*xi2/(r2(i)**2))*(signyy(i)-xi2*signxx(i))
341 . + two*(alpha3(i)/(r1c(i)**2))*signxx(i)
342 normyy = -two*(alpha1(i)*xi1/(r1(i)**2))*(signxx(i)-xi1*signyy(i))
343 . + two*(alpha2(i)/(r2(i)**2))*(signyy(i)-xi2*signxx(i))
344 . + two*(alpha4(i)/(r2c(i)**2))*signyy(i)
345 normxy = (alpha5(i)/(rt(i)**2))*abs(signxy(i))*sign(one,signxy(i))
346c
347 ! 2 - Computation of DFP_DSIG the normal to the non-associated plastic potential
348 !-------------------------------------------------------------------------------
349 ! a) Computation of derivatives
350 normpxx = beta1(i)*(two*signxx(i) + k2*signyy(i)) + beta2(i)*(two*signxx(i) + k4*signyy(i))
351 normpyy = beta1(i)*(two*k1*signyy(i) + k2*signxx(i)) + beta2(i)*(two*k3*signyy(i) + k4*signxx(i))
352 normpxy = (beta1(i)*k5 + beta2(i)*k6)*signxy(i)
353 ! b) Computation of the norm of the derivative
354 normsig = sqrt(normpxx*normpxx + normpyy*normpyy + two*normpxy*normpxy)
355 normsig = max(normsig,em20)
356 ! c) Computation of the normal to the plastic potential
357 normpxx = normpxx/normsig
358 normpyy = normpyy/normsig
359 normpxy = normpxy/normsig
360c
361 ! 3 - Computation of DPHI_DLAMBDA
362 !---------------------------------------------------------
363c
364 ! a) Derivative with respect stress increments tensor DSIG
365 ! --------------------------------------------------------
366 dfdsig2 = normxx * (a11*normpxx + a12*normpyy)
367 . + normyy * (a21*normpxx + a22*normpyy)
368 . + two*normxy * two*normpxy * g12
369c
370 ! b) Derivatives with respect to hardening parameters
371 ! ---------------------------------------------------
372c
373 ! i) Derivatives of PHI with respect to the yield stresses
374 dphi_dr1 = -two*alpha1(i)*(((signxx(i)-xi1*signyy(i))**2)/(r1(i)**3))
375 dphi_dr2 = -two*alpha2(i)*(((signyy(i)-xi2*signxx(i))**2)/(r2(i)**3))
376 dphi_dr1c = -two*alpha3(i)*(signxx(i)**2)/(r1c(i)**3)
377 dphi_dr2c = -two*alpha4(i)*(signyy(i)**2)/(r2c(i)**3)
378 dphi_drt = -two*alpha5(i)*(signxy(i)**2)/(rt(i)**3)
379 ! ii) Derivatives of yield stresses with respect to hardening parameters
380 IF (itab == 0) THEN
381 dr1_dp(i) = (one/(cini1 + s1*pla(i,2))) * (one - (s1*pla(i,2) /(cini1 + s1*pla(i,2))))
382 dr2_dp(i) = (one/(cini2 + s2*pla(i,3))) * (one - (s2*pla(i,3) /(cini2 + s2*pla(i,3))))
383 dr1c_dp(i) = (one/(cini1c + s1c*pla(i,4))) * (one - (s1c*pla(i,4)/(cini1c + s1c*pla(i,4))))
384 dr1c_dp(i) = dr1c_dp(i)/sqrt(one - g1c)
385 dr2c_dp(i) = (one/(cini2c + s2c*pla(i,5))) * (one - (s2c*pla(i,5)/(cini2c + s2c*pla(i,5))))
386 drt_dp(i) = eta2(i)*(one/(cinit + st*pla(i,6))) * (one - (st*pla(i,6) /(cinit + st*pla(i,6))))
387 ENDIF
388 ! iii) Assembling derivatives of PHI with respect to hardening parameter
389 hardp(i) = sqrt(alpha1(i)*dr1_dp(i)**2 + alpha2(i)*dr2_dp(i)**2 + alpha3(i)*dr1c_dp(i)**2
390 . + alpha4(i)*dr2c_dp(i)**2 + two*alpha5(i)*drt_dp(i)**2)
391 dphi_dpla = dphi_dr1*dr1_dp(i)*alpha1(i) + dphi_dr2*dr2_dp(i)*alpha2(i) +
392 . dphi_dr1c*dr1c_dp(i)*alpha3(i) + dphi_dr2c*dr2c_dp(i)*alpha4(i) +
393 . dphi_drt*drt_dp(i)*alpha5(i)
394c
395 ! iv) Derivative of PHI with respect to DLAM ( = -DENOM)
396 dphi_dlam = -dfdsig2 + dphi_dpla
397 dphi_dlam = sign(max(abs(dphi_dlam),em20) ,dphi_dlam)
398c
399 ! 4 - Computation of plastic multiplier and variables update
400 !----------------------------------------------------------
401c
402 ! a) Computation of the plastic multiplier
403 dlam = -phi(i) / dphi_dlam
404c
405 ! b) Plastic strains tensor update
406 dpxx = dlam * normpxx
407 dpyy = dlam * normpyy
408 dpxy = two * dlam * normpxy
409c
410 ! c) Elasto-plastic stresses update
411 signxx(i) = signxx(i) - (a11*dpxx + a12*dpyy)
412 signyy(i) = signyy(i) - (a21*dpxx + a22*dpyy)
413 signxy(i) = signxy(i) - g12*dpxy
414c
415 ! d) Cumulated plastic strain and hardening parameter update
416 ddep = dlam
417 dpla(i) = max(zero, dpla(i) + ddep)
418 dpla2(i) = max(zero, dpla2(i) + alpha1(i)*ddep)
419 dpla3(i) = max(zero, dpla3(i) + alpha2(i)*ddep)
420 dpla4(i) = max(zero, dpla4(i) + alpha3(i)*ddep)
421 dpla5(i) = max(zero, dpla5(i) + alpha4(i)*ddep)
422 dpla6(i) = max(zero, dpla6(i) + alpha5(i)*ddep)
423 pla(i,1) = pla(i,1) + ddep
424 pla(i,2) = pla(i,2) + alpha1(i)*ddep
425 pla(i,3) = pla(i,3) + alpha2(i)*ddep
426 pla(i,4) = pla(i,4) + alpha3(i)*ddep
427 pla(i,5) = pla(i,5) + alpha4(i)*ddep
428 pla(i,6) = pla(i,6) + alpha5(i)*ddep
429c
430 ! e) Update of Alpha coefficient
431 alpha1(i) = zero
432 alpha2(i) = zero
433 alpha3(i) = zero
434 alpha4(i) = zero
435 alpha5(i) = zero
436 IF ((signxx(i) - xi1*signyy(i)) > zero) alpha1(i) = one
437 IF ((signyy(i) - xi2*signxx(i)) > zero) alpha2(i) = one
438 IF (-signxx(i) > zero) alpha3(i) = one
439 IF (-signyy(i) > zero) alpha4(i) = one
440 IF (abs(signxy(i)) > zero) alpha5(i) = one
441c
442 ! f) Yield stresses update
443 IF (itab == 0) THEN
444 ! i) Tensile yield stresses update
445 ! - in direction 1 (MD)
446 r1(i) = sigy1 + (one/(cini1 + s1*pla(i,2)))*pla(i,2)
447 ! - in direction 2 (CD)
448 r2(i) = sigy2 + (one/(cini2 + s2*pla(i,3)))*pla(i,3)
449 ! ii) Compressive yield stresses update
450 ! - in direction 1 (MD)
451 r1c(i) = sigy1c + (one/(cini1c + s1c*pla(i,4)))*pla(i,4)
452 r1c(i) = r1c(i)/sqrt(one - g1c)
453 ! - in direction 2 (CD)
454 r2c(i) = sigy2c + (one/(cini2c + s2c*pla(i,5)))*pla(i,5)
455 ! iii) Shear yield stress update
456 eta1(i) = one + alpha3(i)*alpha4(i)*d1
457 eta2(i) = one + alpha3(i)*alpha4(i)*d2
458 rt(i) = eta1(i)*sigyt + eta2(i)*(one/(cinit + st*pla(i,6)))*pla(i,6)
459c
460 ! g) Yield function value update
461 phi(i) = alpha1(i)*((signxx(i) - xi1*signyy(i))/r1(i))**2 +
462 . alpha2(i)*((signyy(i) - xi2*signxx(i))/r2(i))**2 +
463 . alpha3(i)*((- signxx(i))/r1c(i))**2 +
464 . alpha4(i)*((- signyy(i))/r2c(i))**2 +
465 . alpha5(i)*(abs(signxy(i))/rt(i))**2 - one
466 ENDIF
467c
468 ENDDO
469 ! End of the loop over yielding elements
470c
471 ! If tabulated yield function, update of the yield stress for all element
472 IF (itab > 0) THEN
473 IF (nindx > 0) THEN
474 ! -> Tensile yield stress in direction 1 (MD)
475 xvec(1:nel,1) = pla(1:nel,2)
476 xvec(1:nel,2) = epsd(1:nel,2) * xscale1
477 ipos(1:nel,1) = vartmp(1:nel,1)
478 ipos(1:nel,2) = 1
479 CALL table2d_vinterp_log(table(itable(1)),ismooth,nel,nel,ipos,xvec,r1,dr1_dp,hardr)
480 r1(1:nel) = r1(1:nel) * yscale1
481 dr1_dp(1:nel) = dr1_dp(1:nel) * yscale1
482 vartmp(1:nel,1) = ipos(1:nel,1)
483 ! -> Compressive yield stress in direction 1 (MD)
484 xvec(1:nel,1) = pla(1:nel,3)
485 xvec(1:nel,2) = epsd(1:nel,3) * xscale2
486 ipos(1:nel,1) = vartmp(1:nel,2)
487 ipos(1:nel,2) = 1
488 CALL table2d_vinterp_log(table(itable(2)),ismooth,nel,nel,ipos,xvec,r2,dr2_dp,hardr)
489 r2(1:nel) = r2(1:nel) * yscale2
490 dr2_dp(1:nel) = dr2_dp(1:nel) * yscale2
491 vartmp(1:nel,2) = ipos(1:nel,1)
492 ! -> Positive shear yield stress
493 xvec(1:nel,1) = pla(1:nel,4)
494 xvec(1:nel,2) = epsd(1:nel,4) * xscale1c
495 ipos(1:nel,1) = vartmp(1:nel,3)
496 ipos(1:nel,2) = 1
497 CALL table2d_vinterp_log(table(itable(3)),ismooth,nel,nel,ipos,xvec,r1c,dr1c_dp,hardr)
498 r1c(1:nel) = r1c(1:nel) * yscale1c
499 dr1c_dp(1:nel)= dr1c_dp(1:nel) * yscale1c
500 vartmp(1:nel,3) = ipos(1:nel,1)
501 ! -> Tensile yield stress in direction 2 (CD)
502 xvec(1:nel,1) = pla(1:nel,5)
503 xvec(1:nel,2) = epsd(1:nel,5) * xscale2c
504 ipos(1:nel,1) = vartmp(1:nel,4)
505 ipos(1:nel,2) = 1
506 CALL table2d_vinterp_log(table(itable(4)),ismooth,nel,nel,ipos,xvec,r2c,dr2c_dp,hardr)
507 r2c(1:nel) = r2c(1:nel) * yscale2c
508 dr2c_dp(1:nel)= dr2c_dp(1:nel) * yscale2c
509 vartmp(1:nel,4) = ipos(1:nel,1)
510 ! Compressive yield stress in direction 2 (CD)
511 xvec(1:nel,1) = pla(1:nel,6)
512 xvec(1:nel,2) = epsd(1:nel,6) * xscalet
513 ipos(1:nel,1) = vartmp(1:nel,5)
514 ipos(1:nel,2) = 1
515 CALL table2d_vinterp_log(table(itable(5)),ismooth,nel,nel,ipos,xvec,rt,drt_dp,hardr)
516 rt(1:nel) = rt(1:nel) * yscalet
517 drt_dp(1:nel) = drt_dp(1:nel) * yscalet
518 vartmp(1:nel,5) = ipos(1:nel,1)
519 ! Yield function value update
520 DO i = 1, nel
521 phi(i) = alpha1(i)*((signxx(i) - xi1*signyy(i))/r1(i))**2 +
522 . alpha2(i)*((signyy(i) - xi2*signxx(i))/r2(i))**2 +
523 . alpha3(i)*((- signxx(i))/r1c(i))**2 +
524 . alpha4(i)*((- signyy(i))/r2c(i))**2 +
525 . alpha5(i)*(abs(signxy(i))/rt(i))**2 - one
526 ENDDO
527 ENDIF
528 ENDIF
529 ENDDO
530 ! End of the loop over the iterations
531 !===================================================================
532 ! - END OF PLASTIC CORRECTION WITH NEWTON IMPLICIT ITERATIVE METHOD
533 !===================================================================
534c
535 ! Storing new values
536 DO i=1,nel
537 ! Plastic strain-rate
538 IF (itab == 0) THEN
539 epsd(i,1) = dpla(i) / max(em20,timestep)
540 epsd(i,2) = dpla2(i) / max(em20,timestep)
541 epsd(i,3) = dpla3(i) / max(em20,timestep)
542 epsd(i,4) = dpla4(i) / max(em20,timestep)
543 epsd(i,5) = dpla5(i) / max(em20,timestep)
544 epsd(i,6) = dpla6(i) / max(em20,timestep)
545 ELSE
546 dpdt = dpla(i)/max(em20,timestep)
547 epsd(i,1) = asrate * dpdt + (one - asrate) * epsd(i,1)
548 dpdt = dpla2(i)/max(em20,timestep)
549 epsd(i,2) = asrate * dpdt + (one - asrate) * epsd(i,2)
550 dpdt = dpla3(i)/max(em20,timestep)
551 epsd(i,3) = asrate * dpdt + (one - asrate) * epsd(i,3)
552 dpdt = dpla4(i)/max(em20,timestep)
553 epsd(i,4) = asrate * dpdt + (one - asrate) * epsd(i,4)
554 dpdt = dpla5(i)/max(em20,timestep)
555 epsd(i,5) = asrate * dpdt + (one - asrate) * epsd(i,5)
556 dpdt = dpla6(i)/max(em20,timestep)
557 epsd(i,6) = asrate * dpdt + (one - asrate) * epsd(i,6)
558 ENDIF
559 ! Coefficient for hourglass
560 et(i) = hardp(i) / (hardp(i) + max(young1,young2))
561 ! Computation of the sound speed
562 soundsp(i) = sqrt(max(a11,a12,a21,a22,young3,g12,g23,g31)/ rho(i))
563 ! Storing the yield stress
564 sigy(i) = sqrt(r1(i)**2 + r2(i)**2 + r1c(i)**2 +
565 . r2c(i)**2 + two*rt(i)**2)
566 ENDDO
567c
my_real
#define my_real
Definition cppsort.cpp:32
alpha2
#define alpha2
Definition eval.h:48
interface_table_mod::table2d_vinterp_log
Definition table_mod.F:261
max
#define max(a, b)
Definition macros.h:21
interface_table_mod
Definition table_mod.F:238
table_mod
Definition table_mod.F:112