OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat107c_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 mat107c_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

◆ mat107c_nice()

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