OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat112_xia_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 mat112_xia_newton (nel, ngl, nuparam, nuvar, time, timestep, uparam, uvar, jthe, off, rho0, rho, pla, dpla, epsd, soundsp, epszz, 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

◆ mat112_xia_newton()

subroutine mat112_xia_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) epszz,
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) vartmp,
integer, dimension(numtabl) itable,
type(ttable), dimension(ntable) table )

Definition at line 33 of file mat112_xia_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,rho0,epszz,
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,4),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),NINDX2,INDEX2(NEL),NINDX3,INDEX3(NEL),
86 . NITER,ITER,ITAB,ISMOOTH,IPOS(NEL,2)
87c
88 my_real
89 . young1,young2,young3,nu12,nu21,
90 . a11,a12,a21,a22,
91 . g12,g23,g31,deuxk,e3c,cc,c1,
92 . nu1p,nu2p,nu4p,nu5p,s01,a01,b01,c01,
93 . s02,a02,b02,c02,s03,a03,b03,c03,
94 . s04,a04,b04,c04,s05,a05,b05,c05,
95 . asig,bsig,csig,tau0,atau,btau,n3,n6,
96 . bulk,nu13,nu31,nu23,nu32
97 my_real
98 . normsig,h,dpdt,dlam,ddep,depxx,depyy
99 my_real
100 . normxx,normyy,normxy,normpxx,normpyy,normpxy,
101 . dphi_dlam,dphi,dfdsig2,dphi_dpla,dpxx,dpyy,dpxy,
102 . dphi_dsigy1,dphi_dsigy2,dphi_dsigy3,dphi_dsigy4,
103 . dphi_dsigy5
104 my_real
105 . normzz,dtheta_dlam,dtheta,dtheta_dpla,
106 . dtheta_dsigyo
107 my_real
108 . normyz,normzx,dpsi_dlam,dpsi,dpsi_dpla,
109 . dpyz,dpzx,dpsi_dsigys
110 my_real
111 . xscale1,yscale1,xscale2,yscale2,xscale3,yscale3,xscale4,yscale4,
112 . xscale5,yscale5,xscalec,yscalec,xscales,yscales,xvec(nel,2),asrate
113c
114 my_real, DIMENSION(2) ::
115 . n1,n2,n4,n5
116c
117 my_real, DIMENSION(NEL) ::
118 . khi1,khi2,khi3,khi4,khi5,khi6,sigy1,sigy2,sigy3,sigy4,sigy5,dpla2,dpla3,
119 . dpla4,sigys,sigyo,dsigxx,dsigyy,dsigxy,dsigyz,dsigzx,hardp,phi,phi0,psi,psi0,
120 . theta,theta0,dsigzz,eezz,gam1,gam2,gam3,gam4,gam5,gam6,dsigy1_dp,dsigy2_dp,
121 . dsigy3_dp,dsigy4_dp,dsigy5_dp,dsigyo_dp,dsigys_dp,hardr
122c
123 !=======================================================================
124 ! - INITIALISATION OF COMPUTATION ON TIME STEP
125 !=======================================================================
126 ! Recovering model parameters
127 ! Elastic parameters
128 young1 = uparam(1) ! Young modulus in direction 1 (MD)
129 young2 = uparam(2) ! Young modulus in direction 2 (CD)
130 young3 = uparam(3) ! young modulus in direction 3 (out of plane)
131 nu12 = uparam(4) ! Poisson's ratio in 12
132 nu21 = uparam(5) ! Poisson's ratio in 21
133 a11 = uparam(6) ! Component 11 of orthotropic shell elasticity matrix
134 a12 = uparam(7) ! Component 12 of orthotropic shell elasticity matrix
135 a21 = uparam(8) ! Component 21 of orthotropic shell elasticity matrix
136 a22 = uparam(9) ! Component 22 of orthotropic shell elasticity matrix
137 g12 = uparam(10) ! Shear modulus in 12
138 g23 = uparam(11) ! Shear modulus in 23
139 g31 = uparam(12) ! Shear modulus in 31
140 itab = int(uparam(14)) ! Tabulated yield stress flag
141 deuxk = uparam(15) ! Yield criterion exponent
142 e3c = uparam(16) ! Compressive elasticity parameter
143 cc = uparam(17) ! Exponent compressive elasticity parameter
144 nu1p = uparam(18) ! Tensile plastic Poisson ratio in direction 1 (MD)
145 nu2p = uparam(19) ! Tensile plastic Poisson ratio in direction 2 (CD)
146 nu4p = uparam(20) ! Compressive plastic Poisson ratio in direction 1 (MD)
147 nu5p = uparam(21) ! Compressive plastic Poisson ratio in direction 2 (CD)
148 IF (itab == 0) THEN
149 s01 = uparam(22) ! 1st Tensile plasticity parameter direction 1 (MD)
150 a01 = uparam(23) ! 2nd Tensile plasticity parameter direction 1 (MD)
151 b01 = uparam(24) ! 3rd Tensile plasticity parameter direction 1 (MD)
152 c01 = uparam(25) ! 4th Tensile plasticity parameter direction 1 (MD)
153 s02 = uparam(26) ! 1st Tensile plasticity parameter direction 2 (CD)
154 a02 = uparam(27) ! 2nd Tensile plasticity parameter direction 2 (CD)
155 b02 = uparam(28) ! 3rd Tensile plasticity parameter direction 2 (CD)
156 c02 = uparam(29) ! 4th Tensile plasticity parameter direction 2 (CD)
157 s03 = uparam(30) ! 1st Positive shear plasticity parameter
158 a03 = uparam(31) ! 2nd Positive shear plasticity parameter
159 b03 = uparam(32) ! 3rd Positive shear plasticity parameter
160 c03 = uparam(33) ! 4th Positive shear plasticity parameter
161 s04 = uparam(34) ! 1st Compressive plasticity parameter direction 1 (MD)
162 a04 = uparam(35) ! 2nd Compressive plasticity parameter direction 1 (MD)
163 b04 = uparam(36) ! 3rd Compressive plasticity parameter direction 1 (MD)
164 c04 = uparam(37) ! 4th Compressive plasticity parameter direction 1 (MD)
165 s05 = uparam(38) ! 1st Compressive plasticity parameter direction 2 (CD)
166 a05 = uparam(39) ! 2nd Compressive plasticity parameter direction 2 (CD)
167 b05 = uparam(40) ! 3rd Compressive plasticity parameter direction 2 (CD)
168 c05 = uparam(41) ! 4th Compressive plasticity parameter direction 2 (CD)
169 asig = uparam(42) ! 1st Out-of-plane plasticity parameter
170 bsig = uparam(43) ! 2nd Out-of-plane plasticity parameter
171 csig = uparam(44) ! 3rd Out-of-plane plasticity parameter
172 tau0 = uparam(45) ! 1st Transverse shear plasticity parameter
173 atau = uparam(46) ! 2nd Transverse shear plasticity parameter
174 btau = uparam(47) ! 3rd Transverse shear plasticity parameter
175 ELSE
176 xscale1 = uparam(22)
177 yscale1 = uparam(23)
178 xscale2 = uparam(24)
179 yscale2 = uparam(25)
180 xscale3 = uparam(26)
181 yscale3 = uparam(27)
182 xscale4 = uparam(28)
183 yscale4 = uparam(29)
184 xscale5 = uparam(30)
185 yscale5 = uparam(31)
186 xscalec = uparam(32)
187 yscalec = uparam(33)
188 xscales = uparam(34)
189 yscales = uparam(35)
190 asrate = uparam(36)
191 asrate = (two*pi*asrate*timestep)/(two*pi*asrate*timestep + one)
192 ismooth = int(uparam(37))
193 ENDIF
194c
195 ! Yield planes normal
196 n1(1) = one/(sqrt(one+nu1p**2))
197 n1(2) = -nu1p/(sqrt(one+nu1p**2))
198 n2(1) = -nu2p/(sqrt(one+nu2p**2))
199 n2(2) = one/(sqrt(one+nu2p**2))
200 n3 = one
201 n4(1) = -one/(sqrt(one+nu4p**2))
202 n4(2) = nu4p/(sqrt(one+nu4p**2))
203 n5(1) = nu5p/(sqrt(one+nu5p**2))
204 n5(2) = -one/(sqrt(one+nu5p**2))
205 n6 = -one
206c
207 ! Recovering internal variables
208 DO i=1,nel
209 ! OFF parameter for element deletion
210 IF (off(i) < em01) off(i) = zero
211 IF (off(i) < one) off(i) = off(i)*four_over_5
212 ! User inputs
213 eezz(i) = epszz(i) + pla(i,3) ! Out-of-plane elastic strain
214 ! Standard inputs
215 dpla(i) = zero ! Initialization of the "global" plastic strain increment
216 dpla2(i) = zero ! Initialization of the in-plane plastic strain increment
217 dpla3(i) = zero ! Initialization of the out-of-plane plastic strain increment
218 dpla4(i) = zero ! Initialization of the transverse shear plastic strain increment
219 et(i) = one ! Initialization of hourglass coefficient
220 hardp(i) = zero ! Initialization of hourglass coefficient
221 ENDDO
222c
223 ! Computation of the initial yield stress
224 IF (itab > 0) THEN
225 ! In-plane tabulation with strain-rate
226 xvec(1:nel,1) = pla(1:nel,2)
227 xvec(1:nel,2) = epsd(1:nel,2) * xscale1
228 ! -> Tensile yield stress in direction 1 (MD)
229 ipos(1:nel,1) = vartmp(1:nel,1)
230 ipos(1:nel,2) = 1
231 CALL table2d_vinterp_log(table(itable(1)),ismooth,nel,nel,ipos,xvec,sigy1,dsigy1_dp,hardr)
232 sigy1(1:nel) = sigy1(1:nel) * yscale1
233 dsigy1_dp(1:nel)= dsigy1_dp(1:nel) * yscale1
234 vartmp(1:nel,1) = ipos(1:nel,1)
235 ! -> Tensile yield stress in direction 2 (CD)
236 xvec(1:nel,2) = epsd(1:nel,2) * xscale2
237 ipos(1:nel,1) = vartmp(1:nel,2)
238 ipos(1:nel,2) = 1
239 CALL table2d_vinterp_log(table(itable(2)),ismooth,nel,nel,ipos,xvec,sigy2,dsigy2_dp,hardr)
240 sigy2(1:nel) = sigy2(1:nel) * yscale2
241 dsigy2_dp(1:nel)= dsigy2_dp(1:nel) * yscale2
242 vartmp(1:nel,2) = ipos(1:nel,1)
243 ! -> Positive shear yield stress
244 xvec(1:nel,2) = epsd(1:nel,2) * xscale3
245 ipos(1:nel,1) = vartmp(1:nel,3)
246 ipos(1:nel,2) = 1
247 CALL table2d_vinterp_log(table(itable(3)),ismooth,nel,nel,ipos,xvec,sigy3,dsigy3_dp,hardr)
248 sigy3(1:nel) = sigy3(1:nel) * yscale3
249 dsigy3_dp(1:nel)= dsigy3_dp(1:nel) * yscale3
250 vartmp(1:nel,3) = ipos(1:nel,1)
251 ! -> Compressive yield stress in direction 1 (MD)
252 xvec(1:nel,2) = epsd(1:nel,2) * xscale4
253 ipos(1:nel,1) = vartmp(1:nel,4)
254 ipos(1:nel,2) = 1
255 CALL table2d_vinterp_log(table(itable(4)),ismooth,nel,nel,ipos,xvec,sigy4,dsigy4_dp,hardr)
256 sigy4(1:nel) = sigy4(1:nel) * yscale4
257 dsigy4_dp(1:nel)= dsigy4_dp(1:nel) * yscale4
258 vartmp(1:nel,4) = ipos(1:nel,1)
259 ! -> Compressive yield stress in direction 2 (CD)
260 xvec(1:nel,2) = epsd(1:nel,2) * xscale5
261 ipos(1:nel,1) = vartmp(1:nel,5)
262 ipos(1:nel,2) = 1
263 CALL table2d_vinterp_log(table(itable(5)),ismooth,nel,nel,ipos,xvec,sigy5,dsigy5_dp,hardr)
264 sigy5(1:nel) = sigy5(1:nel) * yscale5
265 dsigy5_dp(1:nel)= dsigy5_dp(1:nel) * yscale5
266 vartmp(1:nel,5) = ipos(1:nel,1)
267 ! -> Out-of-plane tabulation with strain-rate
268 xvec(1:nel,1) = pla(1:nel,3)
269 xvec(1:nel,2) = epsd(1:nel,3) * xscalec
270 ! -> Transverse shear yield stress
271 ipos(1:nel,1) = vartmp(1:nel,6)
272 ipos(1:nel,2) = 1
273 CALL table2d_vinterp_log(table(itable(6)),ismooth,nel,nel,ipos,xvec,sigyo,dsigyo_dp,hardr)
274 sigyo(1:nel) = sigyo(1:nel) * yscalec
275 dsigyo_dp(1:nel)= dsigyo_dp(1:nel) * yscalec
276 vartmp(1:nel,6) = ipos(1:nel,1)
277 ! -> Transverse shear tabulation with strain-rate
278 xvec(1:nel,1) = pla(1:nel,4)
279 xvec(1:nel,2) = epsd(1:nel,4) * xscales
280 ! -> Transverse shear yield stress
281 ipos(1:nel,1) = vartmp(1:nel,7)
282 ipos(1:nel,2) = 1
283 CALL table2d_vinterp_log(table(itable(7)),ismooth,nel,nel,ipos,xvec,sigys,dsigys_dp,hardr)
284 sigys(1:nel) = sigys(1:nel) * yscales
285 dsigys_dp(1:nel)= dsigys_dp(1:nel) * yscales
286 vartmp(1:nel,7) = ipos(1:nel,1)
287 ELSE
288 DO i = 1,nel
289 ! Tensile yield stress in direction 1 (MD)
290 sigy1(i) = s01 + a01*tanh(b01*pla(i,2)) + c01*pla(i,2)
291 ! Tensile yield stress in direction 2 (CD)
292 sigy2(i) = s02 + a02*tanh(b02*pla(i,2)) + c02*pla(i,2)
293 ! Positive shear yield stress
294 sigy3(i) = s03 + a03*tanh(b03*pla(i,2)) + c03*pla(i,2)
295 ! Compressive yield stress in direction 1 (MD)
296 sigy4(i) = s04 + a04*tanh(b04*pla(i,2)) + c04*pla(i,2)
297 ! Compressive yield stress in direction 2 (CD)
298 sigy5(i) = s05 + a05*tanh(b05*pla(i,2)) + c05*pla(i,2)
299 ! Out-of-plane yield stress
300 sigyo(i) = asig + bsig*exp(csig*pla(i,3))
301 ! Transverse shear yield stress
302 sigys(i) = tau0 + (atau - min(zero,sigozz(i))*btau)*pla(i,4)
303 ENDDO
304 ENDIF
305c
306 !========================================================================
307 ! - COMPUTATION OF TRIAL VALUES
308 !========================================================================
309 DO i=1,nel
310c
311 ! Computation of the trial stress tensor
312 ! - in-plane stresses
313 signxx(i) = sigoxx(i) + a11*depsxx(i) + a12*depsyy(i)
314 signyy(i) = sigoyy(i) + a21*depsxx(i) + a22*depsyy(i)
315 signxy(i) = sigoxy(i) + depsxy(i)*g12
316 ! - out-of-plane stress
317 IF (eezz(i) >= zero) THEN
318 signzz(i) = sigozz(i) + young3*depszz(i)
319 ELSE
320 signzz(i) = sigozz(i) + e3c*cc*exp(-cc*eezz(i))*depszz(i)
321 ENDIF
322 ! - transverse shear stresses
323 signyz(i) = sigoyz(i) + depsyz(i)*g23
324 signzx(i) = sigozx(i) + depszx(i)*g31
325c
326 ! Computation of trial khi coefficients
327 khi1(i) = zero
328 khi2(i) = zero
329 khi3(i) = zero
330 khi4(i) = zero
331 khi5(i) = zero
332 khi6(i) = zero
333 IF ((n1(1)*signxx(i)+n1(2)*signyy(i)) > zero) khi1(i) = one
334 IF ((n2(1)*signxx(i)+n2(2)*signyy(i)) > zero) khi2(i) = one
335 IF (n3*signxy(i) > zero) khi3(i) = one
336 IF ((n4(1)*signxx(i)+n4(2)*signyy(i)) > zero) khi4(i) = one
337 IF ((n5(1)*signxx(i)+n5(2)*signyy(i)) > zero) khi5(i) = one
338 IF (n6*signxy(i) > zero) khi6(i) = one
339c
340 ! Computation of the in-plane yield function
341 gam1(i) = (n1(1)*signxx(i)+n1(2)*signyy(i))/sigy1(i)
342 gam2(i) = (n2(1)*signxx(i)+n2(2)*signyy(i))/sigy2(i)
343 gam3(i) = n3*signxy(i)/sigy3(i)
344 gam4(i) = (n4(1)*signxx(i)+n4(2)*signyy(i))/sigy4(i)
345 gam5(i) = (n5(1)*signxx(i)+n5(2)*signyy(i))/sigy5(i)
346 gam6(i) = n6*signxy(i)/sigy3(i)
347 phi(i) = khi1(i)*gam1(i)**deuxk + khi2(i)*gam2(i)**deuxk +
348 . khi3(i)*gam3(i)**deuxk + khi4(i)*gam4(i)**deuxk +
349 . khi5(i)*gam5(i)**deuxk + khi6(i)*gam6(i)**deuxk - one
350c
351 ! Computation of the out-of-plane shear yield function
352 theta(i) = - signzz(i) - sigyo(i)
353c
354 ! Computation of the transverse shear yield function
355 psi(i) = (sqrt(signyz(i)**2 + signzx(i)**2)/(sigys(i))) - one
356c
357 ENDDO
358c
359 !========================================================================
360 ! - CHECKING THE PLASTIC BEHAVIOR
361 !========================================================================
362c
363 ! Checking plastic behavior
364 nindx = 0
365 nindx2 = 0
366 nindx3 = 0
367 DO i=1,nel
368 ! In plane
369 IF (phi(i) > zero .AND. off(i) == one) THEN
370 nindx=nindx+1
371 index(nindx)=i
372 ENDIF
373 ! Out-of-plane
374 IF (theta(i) > zero .AND. off(i) == one) THEN
375 nindx2=nindx2+1
376 index2(nindx2)=i
377 ENDIF
378 ! Transverse shear
379 IF (psi(i) > zero .AND. off(i) == one) THEN
380 nindx3=nindx3+1
381 index3(nindx3)=i
382 ENDIF
383 ENDDO
384c
385 !====================================================================
386 ! - I IN-PLANE PLASTIC CORRECTION WITH CUTTING PLANE (NEWTON RESOLUTION)
387 !====================================================================
388c
389 ! Number of Newton iterations
390 niter = 3
391c
392 ! Loop over the iterations
393 DO iter = 1, niter
394#include "vectorize.inc"
395 ! Loop over yielding elements
396 DO ii=1,nindx
397 i = index(ii)
398c
399 ! Note : in this part, the purpose is to compute for each iteration
400 ! a plastic multiplier allowing to update internal variables to satisfy
401 ! the consistency condition using the cutting plane method
402 ! Its expression at each iteration is : DLAMBDA = - PHI/DPHI_DLAMBDA
403 ! -> PHI : current value of yield function (known)
404 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
405 ! into account of internal variables kinetic :
406 ! hardening parameters
407c
408 ! 1 - Computation of DPHI_DSIG the normal to the yield criterion
409 !-------------------------------------------------------------
410 ! Computation of derivatives to the yield criterion
411 ! - in-plane normal
412 normxx = deuxk*khi1(i)*(gam1(i)**(deuxk-1))*(n1(1)/sigy1(i)) +
413 . deuxk*khi2(i)*(gam2(i)**(deuxk-1))*(n2(1)/sigy2(i)) +
414 . deuxk*khi4(i)*(gam4(i)**(deuxk-1))*(n4(1)/sigy4(i)) +
415 . deuxk*khi5(i)*(gam5(i)**(deuxk-1))*(n5(1)/sigy5(i))
416 normyy = deuxk*khi1(i)*(gam1(i)**(deuxk-1))*(n1(2)/sigy1(i)) +
417 . deuxk*khi2(i)*(gam2(i)**(deuxk-1))*(n2(2)/sigy2(i)) +
418 . deuxk*khi4(i)*(gam4(i)**(deuxk-1))*(n4(2)/sigy4(i)) +
419 . deuxk*khi5(i)*(gam5(i)**(deuxk-1))*(n5(2)/sigy5(i))
420 normxy = deuxk/two*khi3(i)*(gam3(i)**(deuxk-1))*(n3/sigy3(i)) +
421 . deuxk/two*khi6(i)*(gam6(i)**(deuxk-1))*(n6/sigy3(i))
422c
423 ! Normalization of the derivatives
424 ! - in-plane normal
425 normsig = sqrt(normxx**2 + normyy**2 + two*normxy**2)
426 normsig = max(normsig,em20)
427 normpxx = normxx/normsig
428 normpyy = normyy/normsig
429 normpxy = normxy/normsig
430c
431 ! 3 - Computation of DPHI_DLAMBDA
432 !---------------------------------------------------------
433c
434 ! a) Derivative with respect stress increments tensor DSIG
435 ! --------------------------------------------------------
436 dfdsig2 = normxx * (a11*normpxx + a12*normpyy)
437 . + normyy * (a21*normpxx + a22*normpyy)
438 . + two*normxy * two*normpxy * g12
439c
440 ! b) Derivatives with respect to hardening parameters
441 ! ---------------------------------------------------
442c
443 ! i) Derivatives of PHI with respect to the yield stresses
444 dphi_dsigy1 = deuxk*khi1(i)*(gam1(i)**(deuxk-1))*(-(n1(1)*signxx(i)+n1(2)*signyy(i))/(sigy1(i)**2))
445 dphi_dsigy2 = deuxk*khi2(i)*(gam2(i)**(deuxk-1))*(-(n2(1)*signxx(i)+n2(2)*signyy(i))/(sigy2(i)**2))
446 dphi_dsigy3 = deuxk*khi3(i)*(gam3(i)**(deuxk-1))*(-n3*signxy(i)/(sigy3(i)**2)) +
447 . deuxk*khi6(i)*(gam6(i)**(deuxk-1))*(-n6*signxy(i)/(sigy3(i)**2))
448 dphi_dsigy4 = deuxk*khi4(i)*(gam4(i)**(deuxk-1))*(-(n4(1)*signxx(i)+n4(2)*signyy(i))/(sigy4(i)**2))
449 dphi_dsigy5 = deuxk*khi5(i)*(gam5(i)**(deuxk-1))*(-(n5(1)*signxx(i)+n5(2)*signyy(i))/(sigy5(i)**2))
450 ! ii) Derivatives of yield stresses with respect to hardening parameters
451 IF (itab == 0) THEN
452 dsigy1_dp(i) = a01*b01*(one-(tanh(b01*pla(i,2)))**2) + c01
453 dsigy2_dp(i) = a02*b02*(one-(tanh(b02*pla(i,2)))**2) + c02
454 dsigy3_dp(i) = a03*b03*(one-(tanh(b03*pla(i,2)))**2) + c03
455 dsigy4_dp(i) = a04*b04*(one-(tanh(b04*pla(i,2)))**2) + c04
456 dsigy5_dp(i) = a05*b05*(one-(tanh(b05*pla(i,2)))**2) + c05
457 ENDIF
458 ! iii) Assembling derivatives of PHI with respect to hardening parameter
459 hardp(i) = sqrt(khi1(i)*dsigy1_dp(i)**2 + khi2(i)*dsigy2_dp(i)**2 +
460 . khi3(i)*two*dsigy3_dp(i)**2 + khi4(i)*dsigy4_dp(i)**2 +
461 . khi5(i)*dsigy5_dp(i)**2)
462 dphi_dpla = dphi_dsigy1*dsigy1_dp(i) + dphi_dsigy2*dsigy2_dp(i) +
463 . dphi_dsigy3*dsigy3_dp(i) + dphi_dsigy4*dsigy4_dp(i) +
464 . dphi_dsigy5*dsigy5_dp(i)
465c
466 ! iv) Derivative of PHI with respect to DLAM ( = -DENOM)
467 dphi_dlam = -dfdsig2 + dphi_dpla
468 dphi_dlam = sign(max(abs(dphi_dlam),em20) ,dphi_dlam)
469c
470 ! 4 - Computation of plastic multiplier and variables update
471 !----------------------------------------------------------
472c
473 ! a) Computation of the plastic multiplier
474 dlam = -phi(i) / dphi_dlam
475c
476 ! b) Plastic strains tensor update
477 dpxx = dlam * normpxx
478 dpyy = dlam * normpyy
479 dpxy = two * dlam * normpxy
480c
481 ! c) Elasto-plastic stresses update
482 signxx(i) = signxx(i) - (a11*dpxx + a12*dpyy)
483 signyy(i) = signyy(i) - (a21*dpxx + a22*dpyy)
484 signxy(i) = signxy(i) - g12*dpxy
485c
486 ! d) Cumulated plastic strain and hardening parameter update
487 ddep = dlam
488 dpla2(i) = max(zero, dpla2(i) + ddep)
489 pla(i,2) = pla(i,2) + ddep
490c
491 ! e) Yield stresses update
492 IF (itab == 0) THEN
493 ! Tensile yield stress in direction 1 (MD)
494 sigy1(i) = s01 + a01*tanh(b01*pla(i,2)) + c01*pla(i,2)
495 ! Tensile yield stress in direction 2 (CD)
496 sigy2(i) = s02 + a02*tanh(b02*pla(i,2)) + c02*pla(i,2)
497 ! Positive shear yield stress
498 sigy3(i) = s03 + a03*tanh(b03*pla(i,2)) + c03*pla(i,2)
499 ! Compressive yield stress in direction 1 (MD)
500 sigy4(i) = s04 + a04*tanh(b04*pla(i,2)) + c04*pla(i,2)
501 ! Compressive yield stress in direction 2 (CD)
502 sigy5(i) = s05 + a05*tanh(b05*pla(i,2)) + c05*pla(i,2)
503 ENDIF
504c
505 ! i) Update of trial alpha coefficients
506 khi1(i) = zero
507 khi2(i) = zero
508 khi3(i) = zero
509 khi4(i) = zero
510 khi5(i) = zero
511 khi6(i) = zero
512 IF ((n1(1)*signxx(i)+n1(2)*signyy(i)) > zero) khi1(i) = one
513 IF ((n2(1)*signxx(i)+n2(2)*signyy(i)) > zero) khi2(i) = one
514 IF (n3*signxy(i) > zero) khi3(i) = one
515 IF ((n4(1)*signxx(i)+n4(2)*signyy(i)) > zero) khi4(i) = one
516 IF ((n5(1)*signxx(i)+n5(2)*signyy(i)) > zero) khi5(i) = one
517 IF (n6*signxy(i) > zero) khi6(i) = one
518c
519 ! j) Yield function value update
520 IF (itab == 0) THEN
521 gam1(i) = (n1(1)*signxx(i)+n1(2)*signyy(i))/sigy1(i)
522 gam2(i) = (n2(1)*signxx(i)+n2(2)*signyy(i))/sigy2(i)
523 gam3(i) = n3*signxy(i)/sigy3(i)
524 gam4(i) = (n4(1)*signxx(i)+n4(2)*signyy(i))/sigy4(i)
525 gam5(i) = (n5(1)*signxx(i)+n5(2)*signyy(i))/sigy5(i)
526 gam6(i) = n6*signxy(i)/sigy3(i)
527 phi(i) = khi1(i)*gam1(i)**deuxk + khi2(i)*gam2(i)**deuxk +
528 . khi3(i)*gam3(i)**deuxk + khi4(i)*gam4(i)**deuxk +
529 . khi5(i)*gam5(i)**deuxk + khi6(i)*gam6(i)**deuxk - one
530 ENDIF
531c
532 ENDDO
533 ! End of the loop over yielding elements
534c
535 ! If tabulated yield function, update of the yield stress for all element
536 IF (itab > 0) THEN
537 IF (nindx > 0) THEN
538 ! In-plane tabulation with strain-rate
539 xvec(1:nel,1) = pla(1:nel,2)
540 xvec(1:nel,2) = epsd(1:nel,2) * xscale1
541 ! -> Tensile yield stress in direction 1 (MD)
542 ipos(1:nel,1) = vartmp(1:nel,1)
543 ipos(1:nel,2) = 1
544 CALL table2d_vinterp_log(table(itable(1)),ismooth,nel,nel,ipos,xvec,sigy1,dsigy1_dp,hardr)
545 sigy1(1:nel) = sigy1(1:nel) * yscale1
546 dsigy1_dp(1:nel)= dsigy1_dp(1:nel) * yscale1
547 vartmp(1:nel,1) = ipos(1:nel,1)
548 ! -> Tensile yield stress in direction 2 (CD)
549 xvec(1:nel,2) = epsd(1:nel,2) * xscale2
550 ipos(1:nel,1) = vartmp(1:nel,2)
551 ipos(1:nel,2) = 1
552 CALL table2d_vinterp_log(table(itable(2)),ismooth,nel,nel,ipos,xvec,sigy2,dsigy2_dp,hardr)
553 sigy2(1:nel) = sigy2(1:nel) * yscale2
554 dsigy2_dp(1:nel)= dsigy2_dp(1:nel) * yscale2
555 vartmp(1:nel,2) = ipos(1:nel,1)
556 ! -> Positive shear yield stress
557 xvec(1:nel,2) = epsd(1:nel,2) * xscale3
558 ipos(1:nel,1) = vartmp(1:nel,3)
559 ipos(1:nel,2) = 1
560 CALL table2d_vinterp_log(table(itable(3)),ismooth,nel,nel,ipos,xvec,sigy3,dsigy3_dp,hardr)
561 sigy3(1:nel) = sigy3(1:nel) * yscale3
562 dsigy3_dp(1:nel)= dsigy3_dp(1:nel) * yscale3
563 vartmp(1:nel,3) = ipos(1:nel,1)
564 ! -> Compressive yield stress in direction 1 (MD)
565 xvec(1:nel,2) = epsd(1:nel,2) * xscale4
566 ipos(1:nel,1) = vartmp(1:nel,4)
567 ipos(1:nel,2) = 1
568 CALL table2d_vinterp_log(table(itable(4)),ismooth,nel,nel,ipos,xvec,sigy4,dsigy4_dp,hardr)
569 sigy4(1:nel) = sigy4(1:nel) * yscale4
570 dsigy4_dp(1:nel)= dsigy4_dp(1:nel) * yscale4
571 vartmp(1:nel,4) = ipos(1:nel,1)
572 ! -> Compressive yield stress in direction 2 (CD)
573 xvec(1:nel,2) = epsd(1:nel,2) * xscale5
574 ipos(1:nel,1) = vartmp(1:nel,5)
575 ipos(1:nel,2) = 1
576 CALL table2d_vinterp_log(table(itable(5)),ismooth,nel,nel,ipos,xvec,sigy5,dsigy5_dp,hardr)
577 sigy5(1:nel) = sigy5(1:nel) * yscale5
578 dsigy5_dp(1:nel)= dsigy5_dp(1:nel) * yscale5
579 vartmp(1:nel,5) = ipos(1:nel,1)
580 ! Yield function value update
581#include "vectorize.inc"
582 DO ii=1,nindx
583 i = index(ii)
584 gam1(i) = (n1(1)*signxx(i)+n1(2)*signyy(i))/sigy1(i)
585 gam2(i) = (n2(1)*signxx(i)+n2(2)*signyy(i))/sigy2(i)
586 gam3(i) = n3*signxy(i)/sigy3(i)
587 gam4(i) = (n4(1)*signxx(i)+n4(2)*signyy(i))/sigy4(i)
588 gam5(i) = (n5(1)*signxx(i)+n5(2)*signyy(i))/sigy5(i)
589 gam6(i) = n6*signxy(i)/sigy3(i)
590 phi(i) = khi1(i)*gam1(i)**deuxk + khi2(i)*gam2(i)**deuxk +
591 . khi3(i)*gam3(i)**deuxk + khi4(i)*gam4(i)**deuxk +
592 . khi5(i)*gam5(i)**deuxk + khi6(i)*gam6(i)**deuxk - one
593 ENDDO
594 ENDIF
595 ENDIF
596 ENDDO
597 ! End of the loop over the iterations
598 !============================================================================
599 ! - END OF IN-PLANE PLASTIC CORRECTION WITH NEWTON IMPLICIT ITERATIVE METHOD
600 !============================================================================
601c
602 !============================================================================
603 ! - II OUT-OF-PLANE PLASTIC CORRECTION WITH CUTTING PLANE (NEWTON RESOLUTION)
604 !============================================================================
605c
606 ! Loop over the iterations
607 DO iter = 1, niter
608#include "vectorize.inc"
609 ! Loop over yielding elements
610 DO ii=1,nindx2
611 i = index2(ii)
612c
613 ! Note : in this part, the purpose is to compute for each iteration
614 ! a plastic multiplier allowing to update internal variables to satisfy
615 ! the consistency condition using the cutting plane method
616 ! Its expression at each iteration is : DLAMBDA = - THETA/DTHETA_DLAMBDA
617 ! -> THETA : current value of yield function (known)
618 ! -> DTHETA_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
619 ! into account of internal variables kinetic :
620 ! hardening parameters
621c
622 ! 1 - Computation of DTHETA_DSIG the normal to the yield criterion
623 !-------------------------------------------------------------
624 ! Computation of derivatives to the yield criterion
625 ! - in-plane normal
626 normzz = -one
627c
628 ! 3 - Computation of DTHETA_DLAMBDA
629 !---------------------------------------------------------
630c
631 ! a) Derivative with respect stress increments tensor DSIG
632 ! --------------------------------------------------------
633 IF (eezz(i) >= zero) THEN
634 dfdsig2 = normzz * young3 * normzz
635 ELSE
636 dfdsig2 = normzz * e3c*cc*exp(-cc*eezz(i)) * normzz
637 ENDIF
638c
639 ! b) Derivatives with respect to hardening parameters
640 ! ---------------------------------------------------
641c
642 ! i) Derivatives of THETA with respect to the yield stresses
643 dtheta_dsigyo = -one
644 ! ii) Derivatives of yield stresses with respect to hardening parameters
645 IF (itab == 0) dsigyo_dp(i) = csig*bsig*exp(csig*pla(i,3))
646 ! iii) Assembling derivatives of THETA with respect to hardening parameter
647 dtheta_dpla = dtheta_dsigyo*dsigyo_dp(i)
648c
649 ! iv) Derivative of THETA with respect to DLAM ( = -DENOM)
650 dtheta_dlam = -dfdsig2 + dtheta_dpla
651 dtheta_dlam = sign(max(abs(dtheta_dlam),em20) ,dtheta_dlam)
652c
653 ! 4 - Computation of plastic multiplier and variables update
654 !----------------------------------------------------------
655c
656 ! a) Computation of the plastic multiplier
657 dlam = -theta(i) / dtheta_dlam
658c
659 ! c) Elasto-plastic stresses update
660 IF (eezz(i)>=zero) THEN
661 signzz(i) = signzz(i) - young3*dlam*normzz
662 ELSE
663 signzz(i) = signzz(i) - e3c*cc*exp(-cc*eezz(i))*dlam*normzz
664 ENDIF
665c
666 ! d) Cumulated plastic strain and hardening parameter update
667 ddep = dlam
668 dpla3(i) = max(zero, dpla3(i) + ddep)
669 pla(i,3) = pla(i,3) + ddep
670 eezz(i) = epszz(i) + pla(i,3)
671c
672 ! e) Yield stresses update
673 IF (itab == 0) THEN
674 ! Out-of-plane yield stress
675 sigyo(i) = asig + bsig*exp(csig*pla(i,3))
676 ! Yield function value update
677 theta(i) = - signzz(i) - sigyo(i)
678 ENDIF
679c
680 ENDDO
681 ! End of the loop over yielding elements
682c
683 ! If tabulated yield function, update of the yield stress for all element
684 IF (itab > 0) THEN
685 IF (nindx2 > 0) THEN
686 ! -> Transverse shear tabulation with strain-rate
687 xvec(1:nel,1) = pla(1:nel,3)
688 xvec(1:nel,2) = epsd(1:nel,3) * xscalec
689 ! -> Transverse shear yield stress
690 ipos(1:nel,1) = vartmp(1:nel,6)
691 ipos(1:nel,2) = 1
692 CALL table2d_vinterp_log(table(itable(6)),ismooth,nel,nel,ipos,xvec,sigyo,dsigyo_dp,hardr)
693 sigyo(1:nel) = sigyo(1:nel) * yscalec
694 dsigyo_dp(1:nel)= dsigyo_dp(1:nel) * yscalec
695 vartmp(1:nel,6) = ipos(1:nel,1)
696 ! Yield function value update
697#include "vectorize.inc"
698 DO ii=1,nindx2
699 i = index2(ii)
700 theta(i) = - signzz(i) - sigyo(i)
701 ENDDO
702 ENDIF
703 ENDIF
704 ENDDO
705 ! End of the loop over the iterations
706 !===============================================================================
707 ! - END OF OUT-OF-PLANE PLASTIC CORRECTION WITH NEWTON IMPLICIT ITERATIVE METHOD
708 !===============================================================================
709c
710 !================================================================================
711 ! - III TRANSVERSE SHEAR PLASTIC CORRECTION WITH CUTTING PLANE (NEWTON RESOLUTION)
712 !================================================================================
713c
714 ! Number of Newton iterations
715 niter = 3
716c
717 ! Loop over the iterations
718 DO iter = 1, niter
719#include "vectorize.inc"
720 ! Loop over yielding elements
721 DO ii=1,nindx3
722 i = index3(ii)
723c
724 ! Note : in this part, the purpose is to compute for each iteration
725 ! a plastic multiplier allowing to update internal variables to satisfy
726 ! the consistency condition using the cutting plane method
727 ! Its expression at each iteration is : DLAMBDA = - PSI/DPSI_DLAMBDA
728 ! -> PSI : current value of yield function (known)
729 ! -> DPSI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
730 ! into account of internal variables kinetic :
731 ! hardening parameters
732c
733 ! 1 - Computation of DPSI_DSIG the normal to the yield criterion
734 !-------------------------------------------------------------
735 ! Computation of derivatives to the yield criterion
736 ! - transverse shear normal
737 normyz = signyz(i)/max((sqrt(signyz(i)**2 + signzx(i)**2)*sigys(i)),em20)
738 normzx = signzx(i)/max((sqrt(signyz(i)**2 + signzx(i)**2)*sigys(i)),em20)
739c
740 ! 2 - Computation of DPSI_DLAMBDA
741 !---------------------------------------------------------
742c
743 ! a) Derivative with respect stress increments tensor DSIG
744 ! --------------------------------------------------------
745 dfdsig2 = two*normyz * g23 * two*normyz +
746 . two*normzx * g31 * two*normzx
747c
748 ! b) Derivatives with respect to hardening parameters
749 ! ---------------------------------------------------
750c
751 ! i) Derivatives of PSI with respect to the yield stresses
752 dpsi_dsigys = -sqrt(signyz(i)**2 + signzx(i)**2)/(sigys(i)**2)
753 ! ii) derivatives of yield stresses with respect to hardening parameters
754 IF (itab == 0) dsigys_dp(i) = atau - min(zero,sigozz(i))*btau
755 ! iii) Assembling derivatives of PSI with respect to hardening parameter
756 dpsi_dpla = dpsi_dsigys*dsigys_dp(i)
757c
758 ! iv) Derivative of PSI with respect to DLAM ( = -DENOM)
759 dpsi_dlam = -dfdsig2 + dpsi_dpla
760 dpsi_dlam = sign(max(abs(dpsi_dlam),em20),dpsi_dlam)
761c
762 ! 4 - Computation of plastic multiplier and variables update
763 !----------------------------------------------------------
764c
765 ! a) Computation of the plastic multiplier
766 dlam = -psi(i) / dpsi_dlam
767c
768 ! b) Plastic strains tensor update
769 dpyz = two*dlam * normyz
770 dpzx = two*dlam * normzx
771c
772 ! c) Elasto-plastic stresses update
773 signyz(i) = signyz(i) - g23*dpyz
774 signzx(i) = signzx(i) - g31*dpzx
775c
776 ! d) Cumulated plastic strain and hardening parameter update
777 ddep = dlam
778 dpla4(i) = max(zero, dpla4(i) + ddep)
779 pla(i,4) = pla(i,4) + ddep
780c
781 ! e) Yield stresses update
782 IF (itab == 0) THEN
783 ! Transverse shear yield stress
784 sigys(i) = tau0 + (atau - min(zero,sigozz(i))*btau)*pla(i,4)
785 ! f) Yield function value update
786 psi(i) = (sqrt(signyz(i)**2 + signzx(i)**2)/sigys(i)) - one
787 ENDIF
788c
789 ENDDO
790 ! End of the loop over yielding elements
791c
792 ! If tabulated yield function, update of the yield stress for all element
793 IF (itab > 0) THEN
794 IF (nindx3 > 0) THEN
795 ! Transverse shear tabulation with strain-rate
796 xvec(1:nel,1) = pla(1:nel,4)
797 xvec(1:nel,2) = epsd(1:nel,4) * xscales
798 ! Transverse shear yield stress
799 ipos(1:nel,1) = vartmp(1:nel,7)
800 ipos(1:nel,2) = 1
801 CALL table2d_vinterp_log(table(itable(7)),ismooth,nel,nel,ipos,xvec,sigys,dsigys_dp,hardr)
802 sigys(1:nel) = sigys(1:nel) * yscales
803 dsigys_dp(1:nel)= dsigys_dp(1:nel) * yscales
804 vartmp(1:nel,7) = ipos(1:nel,1)
805 ! Yield function value update
806#include "vectorize.inc"
807 DO ii=1,nindx3
808 i = index3(ii)
809 psi(i) = (sqrt(signyz(i)**2 + signzx(i)**2)/sigys(i)) - one
810 ENDDO
811 ENDIF
812 ENDIF
813 ENDDO
814 ! End of the loop over the iterations
815 !===================================================================================
816 ! - END OF TRANSVERSE SHEAR PLASTIC CORRECTION WITH NEWTON IMPLICIT ITERATIVE METHOD
817 !===================================================================================
818c
819 ! Storing new values
820 DO i=1,nel
821 ! Plastic strain
822 pla(i,1) = sqrt(pla(i,2)**2 + pla(i,3)**2 + pla(i,4)**2)
823 dpla(i) = (pla(i,2)/(max(sqrt(pla(i,2)**2 + pla(i,3)**2 + pla(i,4)**2),em20)))*dpla2(i) +
824 . (pla(i,3)/(max(sqrt(pla(i,2)**2 + pla(i,3)**2 + pla(i,4)**2),em20)))*dpla3(i) +
825 . (pla(i,4)/(max(sqrt(pla(i,2)**2 + pla(i,3)**2 + pla(i,4)**2),em20)))*dpla4(i)
826 ! Plastic strain-rate
827 IF (itab == 0) THEN
828 epsd(i,1) = dpla(i) / max(em20,timestep)
829 epsd(i,2) = dpla2(i) / max(em20,timestep)
830 epsd(i,3) = dpla3(i) / max(em20,timestep)
831 epsd(i,4) = dpla4(i) / max(em20,timestep)
832 ELSE
833 dpdt = dpla(i)/max(em20,timestep)
834 epsd(i,1) = asrate * dpdt + (one - asrate) * epsd(i,1)
835 dpdt = dpla2(i)/max(em20,timestep)
836 epsd(i,2) = asrate * dpdt + (one - asrate) * epsd(i,2)
837 dpdt = dpla3(i)/max(em20,timestep)
838 epsd(i,3) = asrate * dpdt + (one - asrate) * epsd(i,3)
839 dpdt = dpla4(i)/max(em20,timestep)
840 epsd(i,4) = asrate * dpdt + (one - asrate) * epsd(i,4)
841 ENDIF
842 ! Coefficient for hourglass and soundspeed
843 IF (eezz(i) >= zero) THEN
844 IF (dpla(i) > zero) THEN
845 et(i) = hardp(i) / (hardp(i) + max(young1,young2,young3))
846 ELSE
847 et(i) = one
848 ENDIF
849 soundsp(i) = sqrt(max(a11,a12,a21,a22,young3,g12,g23,g31)/ rho(i))
850 ELSE
851 IF (dpla(i) > zero) THEN
852 et(i) = hardp(i) / (hardp(i) + max(young1,young2,e3c*cc*exp(-cc*eezz(i))))
853 ELSE
854 et(i) = one
855 ENDIF
856 soundsp(i) = sqrt(max(a11,a12,a21,a22,e3c*cc*exp(-cc*eezz(i)),g12,g23,g31)/ rho(i))
857 ENDIF
858 ! Storing the yield stress
859 sigy(i) = sqrt(khi1(i)*sigy1(i)**2 + khi2(i)*sigy2(i)**2 + khi4(i)*sigy4(i)**2 +
860 . khi5(i)*sigy5(i)**2 + two*khi3(i)*sigy3(i)**2 + sigyo(i)**2
861 . + two*sigys(i)**2)
862 ENDDO
863c
my_real
#define my_real
Definition cppsort.cpp:32
interface_table_mod::table2d_vinterp_log
Definition table_mod.F:261
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
interface_table_mod
Definition table_mod.F:238
table_mod
Definition table_mod.F:112