OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat107_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 mat107_nice (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)

Function/Subroutine Documentation

◆ mat107_nice()

subroutine mat107_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) rho0,
intent(in) rho,
intent(inout) pla,
intent(inout) dpla,
intent(inout) epsd,
intent(out) soundsp,
intent(in) depsxx,
intent(in) depsyy,
intent(in) depszz,
intent(in) depsxy,
intent(in) depsyz,
intent(in) depszx,
intent(in) sigoxx,
intent(in) sigoyy,
intent(in) sigozz,
intent(in) sigoxy,
intent(in) sigoyz,
intent(in) sigozx,
intent(out) signxx,
intent(out) signyy,
intent(out) signzz,
intent(out) signxy,
intent(out) signyz,
intent(out) signzx,
intent(out) sigy,
intent(out) et,
integer nvartmp,
integer numtabl,
integer, dimension(nel,nvartmp), intent(inout) vartmp,
integer, dimension(numtabl) itable,
type(ttable), dimension(ntable) table )

Definition at line 33 of file mat107_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, 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),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,beta1,beta2,
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,phi0,eta1,eta2,
107 . dr1_dp,dr2_dp,dr1c_dp,dr2c_dp,drt_dp,hardr,dpla2,dpla3,dpla4,
108 . dpla5,dpla6
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 ! User inputs
177 phi0(i) = uvar(i,1) ! Previous yield function value
178 ! Standard inputs
179 dpla(i) = zero ! Initialization of the "global" plastic strain increment
180 dpla2(i) = zero ! Initialization of the in-plane plastic strain increment
181 dpla3(i) = zero ! Initialization of the transverse shear plastic strain increment
182 dpla4(i) = zero ! Initialization of the in-plane plastic strain increment
183 dpla5(i) = zero ! Initialization of the transverse shear plastic strain increment
184 dpla6(i) = zero ! initialization of the in-plane plastic strain increment
185 et(i) = one ! Initialization of hourglass coefficient
186 hardp(i) = zero ! Initialization of hourglass coefficient
187 ENDDO
188c
189 ! Computation of the initial yield stress
190 IF (itab > 0) THEN
191 ! -> Tensile yield stress in direction 1 (MD)
192 xvec(1:nel,1) = pla(1:nel,2)
193 xvec(1:nel,2) = epsd(1:nel,2) * xscale1
194 ipos(1:nel,1) = vartmp(1:nel,1)
195 ipos(1:nel,2) = 1
196 CALL table2d_vinterp_log(table(itable(1)),ismooth,nel,nel,ipos,xvec,r1,dr1_dp,hardr)
197 r1(1:nel) = r1(1:nel) * yscale1
198 dr1_dp(1:nel) = dr1_dp(1:nel) * yscale1
199 vartmp(1:nel,1) = ipos(1:nel,1)
200 ! -> Compressive yield stress in direction 1 (MD)
201 xvec(1:nel,1) = pla(1:nel,3)
202 xvec(1:nel,2) = epsd(1:nel,3) * xscale2
203 ipos(1:nel,1) = vartmp(1:nel,2)
204 ipos(1:nel,2) = 1
205 CALL table2d_vinterp_log(table(itable(2)),ismooth,nel,nel,ipos,xvec,r2,dr2_dp,hardr)
206 r2(1:nel) = r2(1:nel) * yscale2
207 dr2_dp(1:nel) = dr2_dp(1:nel) * yscale2
208 vartmp(1:nel,2) = ipos(1:nel,1)
209 ! -> Positive shear yield stress
210 xvec(1:nel,1) = pla(1:nel,4)
211 xvec(1:nel,2) = epsd(1:nel,4) * xscale1c
212 ipos(1:nel,1) = vartmp(1:nel,3)
213 ipos(1:nel,2) = 1
214 CALL table2d_vinterp_log(table(itable(3)),ismooth,nel,nel,ipos,xvec,r1c,dr1c_dp,hardr)
215 r1c(1:nel) = r1c(1:nel) * yscale1c
216 dr1c_dp(1:nel)= dr1c_dp(1:nel) * yscale1c
217 vartmp(1:nel,3) = ipos(1:nel,1)
218 ! -> Tensile yield stress in direction 2 (CD)
219 xvec(1:nel,1) = pla(1:nel,5)
220 xvec(1:nel,2) = epsd(1:nel,5) * xscale2c
221 ipos(1:nel,1) = vartmp(1:nel,4)
222 ipos(1:nel,2) = 1
223 CALL table2d_vinterp_log(table(itable(4)),ismooth,nel,nel,ipos,xvec,r2c,dr2c_dp,hardr)
224 r2c(1:nel) = r2c(1:nel) * yscale2c
225 dr2c_dp(1:nel)= dr2c_dp(1:nel) * yscale2c
226 vartmp(1:nel,4) = ipos(1:nel,1)
227 ! Compressive yield stress in direction 2 (CD)
228 xvec(1:nel,1) = pla(1:nel,6)
229 xvec(1:nel,2) = epsd(1:nel,6) * xscalet
230 ipos(1:nel,1) = vartmp(1:nel,5)
231 ipos(1:nel,2) = 1
232 CALL table2d_vinterp_log(table(itable(5)),ismooth,nel,nel,ipos,xvec,rt,drt_dp,hardr)
233 rt(1:nel) = rt(1:nel) * yscalet
234 drt_dp(1:nel) = drt_dp(1:nel) * yscalet
235 vartmp(1:nel,5) = ipos(1:nel,1)
236 ENDIF
237c
238 !========================================================================
239 ! - COMPUTATION OF TRIAL VALUES
240 !========================================================================
241 DO i=1,nel
242c
243 ! Computation of the trial stress tensor
244 signxx(i) = sigoxx(i) + a11*depsxx(i) + a12*depsyy(i)
245 signyy(i) = sigoyy(i) + a21*depsxx(i) + a22*depsyy(i)
246 signzz(i) = sigozz(i) + young3*depszz(i)
247 signxy(i) = sigoxy(i) + depsxy(i)*g12
248 signyz(i) = sigoyz(i) + depsyz(i)*g23
249 signzx(i) = sigozx(i) + depszx(i)*g31
250c
251 ! Computation of trial alpha coefficients
252 alpha1(i) = zero
253 alpha2(i) = zero
254 alpha3(i) = zero
255 alpha4(i) = zero
256 alpha5(i) = zero
257 IF ((signxx(i) - xi1*signyy(i)) > zero) alpha1(i) = one
258 IF ((signyy(i) - xi2*signxx(i)) > zero) alpha2(i) = one
259 IF (-signxx(i) > zero) alpha3(i) = one
260 IF (-signyy(i) > zero) alpha4(i) = one
261 IF (abs(signxy(i)) > zero) alpha5(i) = one
262c
263 IF (itab == 0) THEN
264 ! Computation of yield stresses
265 ! Tensile yield stresses
266 ! - in direction 1 (MD)
267 r1(i) = sigy1 + (one/(cini1 + s1*pla(i,2)))*pla(i,2)
268 ! - in direction 2 (CD)
269 r2(i) = sigy2 + (one/(cini2 + s2*pla(i,3)))*pla(i,3)
270 ! Compressive yield stresses
271 ! - in direction 1 (MD) (including correction)
272 r1c(i) = sigy1c + (one/(cini1c + s1c*pla(i,4)))*pla(i,4)
273 r1c(i) = r1c(i)/sqrt(one - g1c)
274 ! - in direction 2 (CD)
275 r2c(i) = sigy2c + (one/(cini2c + s2c*pla(i,5)))*pla(i,5)
276 ! Shear yield stress
277 eta1(i) = one + alpha3(i)*alpha4(i)*d1
278 eta2(i) = one + alpha3(i)*alpha4(i)*d2
279 rt(i) = eta1(i)*sigyt + eta2(i)*(one/(cinit + st*pla(i,6)))*pla(i,6)
280 ENDIF
281c
282 ! Computation of the yield function
283 phi(i) = alpha1(i)*((signxx(i) - xi1*signyy(i))/r1(i))**2 +
284 . alpha2(i)*((signyy(i) - xi2*signxx(i))/r2(i))**2 +
285 . alpha3(i)*((- signxx(i))/r1c(i))**2 +
286 . alpha4(i)*((- signyy(i))/r2c(i))**2 +
287 . alpha5(i)*(abs(signxy(i))/rt(i))**2 - one
288 ENDDO
289c
290 !========================================================================
291 ! - CHECKING THE PLASTIC BEHAVIOR
292 !========================================================================
293c
294 ! Checking plastic behavior for all elements
295 nindx = 0
296 DO i=1,nel
297 IF (phi(i) > zero .AND. off(i) == one) THEN
298 nindx=nindx+1
299 index(nindx)=i
300 ENDIF
301 ENDDO
302c
303 !====================================================================
304 ! - PLASTIC CORRECTION WITH NICE - EXPLICIT 1 METHOD
305 !====================================================================
306#include "vectorize.inc"
307 DO ii=1,nindx
308c
309 ! Number of the element with plastic behaviour
310 i = index(ii)
311c
312 ! Note : in this part, the purpose is to compute for each iteration
313 ! a plastic multiplier allowing to update internal variables to satisfy
314 ! the consistency condition.
315 ! Its expression is : DLAMBDA = - (PHI + DPHI) / DPHI_DLAMBDA
316! -> PHI : current value of yield function (known)
317 ! -> DPHI : yield function prediction (to compute)
318 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
319 ! into account of internal variables kinetic :
320 ! plasticity, temperature and damage (to compute)
321c
322! 1 - Computation of DPHI_DSIG the normal to the yield criterion
323 !-------------------------------------------------------------
324c
325 normxx = two*(alpha1(i)/(r1(i)**2))*(sigoxx(i)-xi1*sigoyy(i))
326 . - two*(alpha2(i)*xi2/(r2(i)**2))*(sigoyy(i)-xi2*sigoxx(i))
327 . + two*(alpha3(i)/(r1c(i)**2))*sigoxx(i)
328 normyy = -two*(alpha1(i)*xi1/(r1(i)**2))*(sigoxx(i)-xi1*sigoyy(i))
329 . + two*(alpha2(i)/(r2(i)**2))*(sigoyy(i)-xi2*sigoxx(i))
330 . + two*(alpha4(i)/(r2c(i)**2))*sigoyy(i)
331 normxy = (alpha5(i)/(rt(i)**2))*abs(sigoxy(i))*sign(one,sigoxy(i))
332c
333! 2 - Computation of DFP_DSIG the normal to the non-associated plastic potential
334 !-------------------------------------------------------------------------------
335 ! a) Computation of BETA switching coefficient
336 normsig = sqrt(sigoxx(i)*sigoxx(i) + sigoyy(i)*sigoyy(i) + two*sigoxy(i)*sigoxy(i))
337 normsig = max(normsig,em20)
338 beta1 = zero
339 beta2 = zero
340 IF ((sigoxx(i)/normsig > em01).OR.(sigoyy(i)/normsig > em01)) beta1 = one
341 IF ((sigoxx(i)/normsig < -em01).OR.(sigoyy(i)/normsig < -em01)) beta2 = one
342 IF (((abs(sigoxx(i))/normsig)<em01).AND.((abs(sigoyy(i))/normsig)<em01)) THEN
343 beta1 = one
344 beta2 = one
345 ENDIF
346 ! b) Computation of derivatives
347 normpxx = beta1*(two*sigoxx(i) + k2*sigoyy(i)) + beta2*(two*sigoxx(i) + k4*sigoyy(i))
348 normpyy = beta1*(two*k1*sigoyy(i) + k2*sigoxx(i)) + beta2*(two*k3*sigoyy(i) + k4*sigoxx(i))
349 normpxy = (beta1*k5 + beta2*k6)*sigoxy(i)
350 ! c) Computation of the norm of the derivative
351 normsig = sqrt(normpxx*normpxx + normpyy*normpyy + two*normpxy*normpxy)
352 normsig = max(normsig,em20)
353 ! d) Computation of the normal to the plastic potential
354 normpxx = normpxx/normsig
355 normpyy = normpyy/normsig
356 normpxy = normpxy/normsig
357c
358 ! 3 - Computation of DPHI_DLAMBDA
359 !---------------------------------------------------------
360c
361 ! a) Derivative with respect stress increments tensor DSIG
362 ! --------------------------------------------------------
363 dfdsig2 = normxx * (a11*normpxx + a12*normpyy)
364 . + normyy * (a21*normpxx + a22*normpyy)
365 . + two*normxy * two*normpxy * g12
366c
367 ! b) Derivatives with respect to hardening parameters
368 ! ---------------------------------------------------
369c
370 ! i) Derivatives of PHI with respect to the yield stresses
371 dphi_dr1 = -two*alpha1(i)*(((sigoxx(i)-xi1*sigoyy(i))**2)/(r1(i)**3))
372 dphi_dr2 = -two*alpha2(i)*(((sigoyy(i)-xi2*sigoxx(i))**2)/(r2(i)**3))
373 dphi_dr1c = -two*alpha3(i)*(sigoxx(i)**2)/(r1c(i)**3)
374 dphi_dr2c = -two*alpha4(i)*(sigoyy(i)**2)/(r2c(i)**3)
375 dphi_drt = -two*alpha5(i)*(sigoxy(i)**2)/(rt(i)**3)
376 ! ii) Derivatives of yield stresses with respect to hardening parameters
377 IF (itab == 0) THEN
378 dr1_dp(i) = (one/(cini1 + s1*pla(i,2))) * (one - (s1*pla(i,2) /(cini1 + s1*pla(i,2))))
379 dr2_dp(i) = (one/(cini2 + s2*pla(i,3))) * (one - (s2*pla(i,3) /(cini2 + s2*pla(i,3))))
380 dr1c_dp(i) = (one/(cini1c + s1c*pla(i,4))) * (one - (s1c*pla(i,4)/(cini1c + s1c*pla(i,4))))
381 dr1c_dp(i) = dr1c_dp(i)/sqrt(one - g1c)
382 dr2c_dp(i) = (one/(cini2c + s2c*pla(i,5))) * (one - (s2c*pla(i,5)/(cini2c + s2c*pla(i,5))))
383 drt_dp(i) = eta2(i)*(one/(cinit + st*pla(i,6))) * (one - (st*pla(i,6) /(cinit + st*pla(i,6))))
384 ENDIF
385 ! iii) Assembling derivatives of PHI with respect to hardening parameter
386 hardp(i) = sqrt(alpha1(i)*dr1_dp(i)**2 + alpha2(i)*dr2_dp(i)**2 + alpha3(i)*dr1c_dp(i)**2
387 . + alpha4(i)*dr2c_dp(i)**2 + two*alpha5(i)*drt_dp(i)**2)
388 dphi_dpla = dphi_dr1*dr1_dp(i)*alpha1(i) + dphi_dr2*dr2_dp(i)*alpha2(i) +
389 . dphi_dr1c*dr1c_dp(i)*alpha3(i) + dphi_dr2c*dr2c_dp(i)*alpha4(i) +
390 . dphi_drt*drt_dp(i)*alpha5(i)
391c
392 ! iv) Derivative of PHI with respect to DLAM ( = -DENOM)
393 dphi_dlam = -dfdsig2 + dphi_dpla
394 dphi_dlam = sign(max(abs(dphi_dlam),em20) ,dphi_dlam)
395c
396 ! 4 - Computation of plastic multiplier and variables update
397 !----------------------------------------------------------
398c
399 ! a) Restoring previous value of the yield function
400 phi(i) = phi0(i)
401c
402 ! b) Computation of the trial stress increment
403 dsigxx(i) = a11*depsxx(i) + a12*depsyy(i)
404 dsigyy(i) = a21*depsxx(i) + a22*depsyy(i)
405 dsigxy(i) = depsxy(i)*g12
406c
407 ! c) Computation of yield surface trial increment DPHI
408 dphi = normxx * dsigxx(i)
409 . + normyy * dsigyy(i)
410 . + two * normxy * dsigxy(i)
411c
412 ! d) computation of the plastic multiplier
413 dlam = -(phi(i) + dphi) / dphi_dlam
414 dlam = max(dlam, zero)
415c
416 ! e) Plastic strains tensor update
417 dpxx = dlam * normpxx
418 dpyy = dlam * normpyy
419 dpxy = two * dlam * normpxy
420c
421 ! f) Elasto-plastic stresses update
422 signxx(i) = sigoxx(i) + dsigxx(i) - (a11*dpxx + a12*dpyy)
423 signyy(i) = sigoyy(i) + dsigyy(i) - (a21*dpxx + a22*dpyy)
424 signxy(i) = sigoxy(i) + dsigxy(i) - g12*dpxy
425c
426 ! g) Cumulated plastic strain and hardening parameter update
427 ddep = dlam
428 dpla(i) = max(zero, dpla(i) + ddep)
429 dpla2(i) = max(zero, dpla2(i) + alpha1(i)*ddep)
430 dpla3(i) = max(zero, dpla3(i) + alpha2(i)*ddep)
431 dpla4(i) = max(zero, dpla4(i) + alpha3(i)*ddep)
432 dpla5(i) = max(zero, dpla5(i) + alpha4(i)*ddep)
433 dpla6(i) = max(zero, dpla6(i) + alpha5(i)*ddep)
434 pla(i,1) = pla(i,1) + ddep
435 pla(i,2) = pla(i,2) + alpha1(i)*ddep
436 pla(i,3) = pla(i,3) + alpha2(i)*ddep
437 pla(i,4) = pla(i,4) + alpha3(i)*ddep
438 pla(i,5) = pla(i,5) + alpha4(i)*ddep
439 pla(i,6) = pla(i,6) + alpha5(i)*ddep
440c
441 ! h) Update of Alpha coefficient
442 alpha1(i) = zero
443 alpha2(i) = zero
444 alpha3(i) = zero
445 alpha4(i) = zero
446 alpha5(i) = zero
447 IF ((signxx(i) - xi1*signyy(i)) > zero) alpha1(i) = one
448 IF ((signyy(i) - xi2*signxx(i)) > zero) alpha2(i) = one
449 IF (-signxx(i) > zero) alpha3(i) = one
450 IF (-signyy(i) > zero) alpha4(i) = one
451 IF (abs(signxy(i)) > zero) alpha5(i) = one
452c
453 ! i) Yield stresses update
454 IF (itab == 0) THEN
455 ! i) Tensile yield stresses update
456 ! - in direction 1 (MD)
457 r1(i) = sigy1 + (one/(cini1 + s1*pla(i,2)))*pla(i,2)
458 ! - in direction 2 (CD)
459 r2(i) = sigy2 + (one/(cini2 + s2*pla(i,3)))*pla(i,3)
460 ! ii) Compressive yield stresses update
461 ! - in direction 1 (MD)
462 r1c(i) = sigy1c + (one/(cini1c + s1c*pla(i,4)))*pla(i,4)
463 r1c(i) = r1c(i)/sqrt(one - g1c)
464 ! - in direction 2 (CD)
465 r2c(i) = sigy2c + (one/(cini2c + s2c*pla(i,5)))*pla(i,5)
466 ! iii) Shear yield stress update
467 eta1(i) = one + alpha3(i)*alpha4(i)*d1
468 eta2(i) = one + alpha3(i)*alpha4(i)*d2
469 rt(i) = eta1(i)*sigyt + eta2(i)*(one/(cinit + st*pla(i,6)))*pla(i,6)
470c
471 ! i) Yield function value update
472 phi(i) = alpha1(i)*((signxx(i) - xi1*signyy(i))/r1(i))**2 +
473 . alpha2(i)*((signyy(i) - xi2*signxx(i))/r2(i))**2 +
474 . alpha3(i)*((- signxx(i))/r1c(i))**2 +
475 . alpha4(i)*((- signyy(i))/r2c(i))**2 +
476 . alpha5(i)*(abs(signxy(i))/rt(i))**2 - one
477 ENDIF
478c
479 ENDDO
480 !===================================================================
481 ! - END OF PLASTIC CORRECTION WITH NICE - EXPLICIT 1 METHOD
482 !===================================================================
483c
484 ! If tabulated yield function, update of the yield stress for all element
485 IF (itab > 0) THEN
486 IF (nindx > 0) THEN
487 ! -> Tensile yield stress in direction 1 (MD)
488 xvec(1:nel,1) = pla(1:nel,2)
489 xvec(1:nel,2) = epsd(1:nel,2) * xscale1
490 ipos(1:nel,1) = vartmp(1:nel,1)
491 ipos(1:nel,2) = 1
492 CALL table2d_vinterp_log(table(itable(1)),ismooth,nel,nel,ipos,xvec,r1,dr1_dp,hardr)
493 r1(1:nel) = r1(1:nel) * yscale1
494 dr1_dp(1:nel) = dr1_dp(1:nel) * yscale1
495 vartmp(1:nel,1) = ipos(1:nel,1)
496 ! -> Compressive yield stress in direction 1 (MD)
497 xvec(1:nel,1) = pla(1:nel,3)
498 xvec(1:nel,2) = epsd(1:nel,3) * xscale2
499 ipos(1:nel,1) = vartmp(1:nel,2)
500 ipos(1:nel,2) = 1
501 CALL table2d_vinterp_log(table(itable(2)),ismooth,nel,nel,ipos,xvec,r2,dr2_dp,hardr)
502 r2(1:nel) = r2(1:nel) * yscale2
503 dr2_dp(1:nel) = dr2_dp(1:nel) * yscale2
504 vartmp(1:nel,2) = ipos(1:nel,1)
505 ! -> Positive shear yield stress
506 xvec(1:nel,1) = pla(1:nel,4)
507 xvec(1:nel,2) = epsd(1:nel,4) * xscale1c
508 ipos(1:nel,1) = vartmp(1:nel,3)
509 ipos(1:nel,2) = 1
510 CALL table2d_vinterp_log(table(itable(3)),ismooth,nel,nel,ipos,xvec,r1c,dr1c_dp,hardr)
511 r1c(1:nel) = r1c(1:nel) * yscale1c
512 dr1c_dp(1:nel)= dr1c_dp(1:nel) * yscale1c
513 vartmp(1:nel,3) = ipos(1:nel,1)
514 ! -> Tensile yield stress in direction 2 (CD)
515 xvec(1:nel,1) = pla(1:nel,5)
516 xvec(1:nel,2) = epsd(1:nel,5) * xscale2c
517 ipos(1:nel,1) = vartmp(1:nel,4)
518 ipos(1:nel,2) = 1
519 CALL table2d_vinterp_log(table(itable(4)),ismooth,nel,nel,ipos,xvec,r2c,dr2c_dp,hardr)
520 r2c(1:nel) = r2c(1:nel) * yscale2c
521 dr2c_dp(1:nel)= dr2c_dp(1:nel) * yscale2c
522 vartmp(1:nel,4) = ipos(1:nel,1)
523 ! compressive yield stress in direction 2 (cd)
524 xvec(1:nel,1) = pla(1:nel,6)
525 xvec(1:nel,2) = epsd(1:nel,6) * xscalet
526 ipos(1:nel,1) = vartmp(1:nel,5)
527 ipos(1:nel,2) = 1
528 CALL table2d_vinterp_log(table(itable(5)),ismooth,nel,nel,ipos,xvec,rt,drt_dp,hardr)
529 rt(1:nel) = rt(1:nel) * yscalet
530 drt_dp(1:nel) = drt_dp(1:nel) * yscalet
531 vartmp(1:nel,5) = ipos(1:nel,1)
532 ! Yield function value update
533 DO i = 1, nel
534 phi(i) = alpha1(i)*((signxx(i) - xi1*signyy(i))/r1(i))**2 +
535 . alpha2(i)*((signyy(i) - xi2*signxx(i))/r2(i))**2 +
536 . alpha3(i)*((- signxx(i))/r1c(i))**2 +
537 . alpha4(i)*((- signyy(i))/r2c(i))**2 +
538 . alpha5(i)*(abs(signxy(i))/rt(i))**2 - one
539 ENDDO
540 ENDIF
541 ENDIF
542c
543 ! Storing new values
544 DO i=1,nel
545 ! USR Outputs
546 uvar(i,1) = phi(i) ! Yield function value
547 ! Plastic strain-rate
548 IF (itab == 0) THEN
549 epsd(i,1) = dpla(i) / max(em20,timestep)
550 epsd(i,2) = dpla2(i) / max(em20,timestep)
551 epsd(i,3) = dpla3(i) / max(em20,timestep)
552 epsd(i,4) = dpla4(i) / max(em20,timestep)
553 epsd(i,5) = dpla5(i) / max(em20,timestep)
554 epsd(i,6) = dpla6(i) / max(em20,timestep)
555 ELSE
556 dpdt = dpla(i)/max(em20,timestep)
557 epsd(i,1) = asrate * dpdt + (one - asrate) * epsd(i,1)
558 dpdt = dpla2(i)/max(em20,timestep)
559 epsd(i,2) = asrate * dpdt + (one - asrate) * epsd(i,2)
560 dpdt = dpla3(i)/max(em20,timestep)
561 epsd(i,3) = asrate * dpdt + (one - asrate) * epsd(i,3)
562 dpdt = dpla4(i)/max(em20,timestep)
563 epsd(i,4) = asrate * dpdt + (one - asrate) * epsd(i,4)
564 dpdt = dpla5(i)/max(em20,timestep)
565 epsd(i,5) = asrate * dpdt + (one - asrate) * epsd(i,5)
566 dpdt = dpla6(i)/max(em20,timestep)
567 epsd(i,6) = asrate * dpdt + (one - asrate) * epsd(i,6)
568 ENDIF
569 ! Coefficient for hourglass
570 et(i) = hardp(i) / (hardp(i) + max(young1,young2))
571 ! Computation of the sound speed
572 soundsp(i) = sqrt(max(a11,a12,a21,a22,young3,g12,g23,g31)/ rho(i))
573 ! Storing the yield stress
574 sigy(i) = sqrt(r1(i)**2 + r2(i)**2 + r1c(i)**2 +
575 . r2c(i)**2 + two*rt(i)**2)
576 ENDDO
577c
my_real
#define my_real
Definition cppsort.cpp:32
the
end diagonal values have been computed in the(sparse) matrix id.SOL
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
interface_table_mod
Definition table_mod.F:238
table_mod
Definition table_mod.F:112