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