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