OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat112c_xia_nice.F File Reference
#include "implicit_f.inc"
#include "com04_c.inc"
#include "vectorize.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine mat112c_xia_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)

Function/Subroutine Documentation

◆ mat112c_xia_nice()

subroutine mat112c_xia_nice ( integer nel,
integer, dimension(nel), intent(in) ngl,
integer nuparam,
integer nuvar,
time,
timestep,
intent(in) uparam,
intent(inout) uvar,
integer jthe,
intent(inout) off,
intent(in) rho,
intent(inout) pla,
intent(inout) dpla,
intent(inout) epsd,
intent(out) soundsp,
intent(in) shf,
intent(in) depsxx,
intent(in) depsyy,
intent(in) depsxy,
intent(in) depsyz,
intent(in) depszx,
intent(in) sigoxx,
intent(in) sigoyy,
intent(in) sigoxy,
intent(in) sigoyz,
intent(in) sigozx,
intent(out) signxx,
intent(out) signyy,
intent(out) signxy,
intent(out) signyz,
intent(out) signzx,
intent(out) sigy,
intent(out) et,
integer numtabl,
integer, dimension(numtabl) itable,
type(ttable), dimension(ntable) table,
integer nvartmp,
integer, dimension(nel,nvartmp) vartmp )

Definition at line 33 of file mat112c_xia_nice.F.

42 !=======================================================================
43 ! Modules
44 !=======================================================================
45 USE table_mod
46 USE interface_table_mod
47 !=======================================================================
48 ! Implicit types
49 !=======================================================================
50#include "implicit_f.inc"
51 !=======================================================================
52 ! Common
53 !=======================================================================
54#include "com04_c.inc"
55 !=======================================================================
56 ! Dummy arguments
57 !=======================================================================
58 INTEGER NEL,NUPARAM,NUVAR,JTHE,NUMTABL,ITABLE(NUMTABL),NVARTMP
59 INTEGER, DIMENSION(NEL), INTENT(IN) :: NGL
60 my_real
61 . time,timestep
62 INTEGER :: VARTMP(NEL,NVARTMP)
63 my_real,DIMENSION(NUPARAM), INTENT(IN) ::
64 . uparam
65 my_real,DIMENSION(NEL), INTENT(IN) ::
66 . rho,shf,
67 . depsxx,depsyy,depsxy,depsyz,depszx,
68 . sigoxx,sigoyy,sigoxy,sigoyz,sigozx
69c
70 my_real ,DIMENSION(NEL), INTENT(OUT) ::
71 . soundsp,sigy,et,
72 . signxx,signyy,signxy,signyz,signzx
73c
74 my_real ,DIMENSION(NEL), INTENT(INOUT) ::
75 . dpla,off
76 my_real ,DIMENSION(NEL,4),INTENT(INOUT) ::
77 . pla,epsd
78 my_real ,DIMENSION(NEL,NUVAR), INTENT(INOUT) ::
79 . uvar
80c
81 TYPE(TTABLE), DIMENSION(NTABLE) :: TABLE
82 !=======================================================================
83 ! Local Variables
84 !=======================================================================
85 INTEGER I,K,II,NINDX,INDEX(NEL),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
my_real
#define my_real
Definition cppsort.cpp:32
interface_table_mod::table2d_vinterp_log
Definition table_mod.F:261
max
#define max(a, b)
Definition macros.h:21
interface_table_mod
Definition table_mod.F:238
table_mod
Definition table_mod.F:112