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