OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat112c_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 mat112c_xia_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

◆ mat112c_xia_newton()

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