OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat112c_xia_newton.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| mat112c_xia_newton ../engine/source/materials/mat/mat112/mat112c_xia_newton.F
25!||--- called by ------------------------------------------------------
26!|| sigeps112c ../engine/source/materials/mat/mat112/sigeps112c.F
27!||--- calls -----------------------------------------------------
28!|| table2d_vinterp_log ../engine/source/tools/curve/table2d_vinterp_log.F
29!||--- uses -----------------------------------------------------
30!|| interface_table_mod ../engine/share/modules/table_mod.F
31!|| table_mod ../engine/share/modules/table_mod.F
32!||====================================================================
33 SUBROUTINE mat112c_xia_newton(
34 1 NEL ,NGL ,NUPARAM ,NUVAR ,TIME ,TIMESTEP,
35 2 UPARAM ,UVAR ,JTHE ,OFF ,RHO ,
36 3 PLA ,DPLA ,EPSD ,SOUNDSP ,SHF ,
37 4 DEPSXX ,DEPSYY ,DEPSXY ,DEPSYZ ,DEPSZX ,
38 5 SIGOXX ,SIGOYY ,SIGOXY ,SIGOYZ ,SIGOZX ,
39 6 SIGNXX ,SIGNYY ,SIGNXY ,SIGNYZ ,SIGNZX ,
40 7 SIGY ,ET ,
41 8 NUMTABL ,ITABLE ,TABLE ,NVARTMP ,VARTMP )
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
698 END
interface_table_mod::table2d_vinterp_log
Definition table_mod.F:261
max
#define max(a, b)
Definition macros.h:21
mat112c_xia_newton
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)
Definition mat112c_xia_newton.F:42
interface_table_mod
Definition table_mod.F:238
table_mod
Definition table_mod.F:112
table_mod::ttable
Definition table_mod.F:125