OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat107c_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 mat107c_newton (nel, ngl, nuparam, nuvar, time, timestep, uparam, uvar, jthe, off, rho, pla, dpla, epsd, soundsp, shf, depsxx, depsyy, depsxy, depsyz, depszx, sigoxx, sigoyy, sigoxy, sigoyz, sigozx, signxx, signyy, signxy, signyz, signzx, sigy, et, numtabl, itable, table, nvartmp, vartmp)

Function/Subroutine Documentation

◆ mat107c_newton()

subroutine mat107c_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) rho,
intent(inout) pla,
intent(inout) dpla,
intent(inout) epsd,
intent(out) soundsp,
intent(in) shf,
intent(in) depsxx,
intent(in) depsyy,
intent(in) depsxy,
intent(in) depsyz,
intent(in) depszx,
intent(in) sigoxx,
intent(in) sigoyy,
intent(in) sigoxy,
intent(in) sigoyz,
intent(in) sigozx,
intent(out) signxx,
intent(out) signyy,
intent(out) signxy,
intent(out) signyz,
intent(out) signzx,
intent(out) sigy,
intent(out) et,
integer numtabl,
integer, dimension(numtabl) itable,
type(ttable), dimension(ntable) table,
integer nvartmp,
integer, dimension(nel,nvartmp) vartmp )

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