OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
mat104c_nodam_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!|| mat104c_nodam_nice ../engine/source/materials/mat/mat104/mat104c_nodam_nice.F
25!||--- called by ------------------------------------------------------
26!|| sigeps104c ../engine/source/materials/mat/mat104/sigeps104c.F
27!||====================================================================
28 SUBROUTINE mat104c_nodam_nice(
29 1 NEL ,NGL ,NUPARAM ,NUVAR ,
30 2 TIME ,TIMESTEP,UPARAM ,UVAR ,JTHE ,OFF ,
31 3 GS ,RHO ,PLA ,DPLA ,EPSD ,SOUNDSP ,
32 4 DEPSXX ,DEPSYY ,DEPSXY ,DEPSYZ ,DEPSZX ,
33 5 SIGOXX ,SIGOYY ,SIGOXY ,SIGOYZ ,SIGOZX ,
34 6 SIGNXX ,SIGNYY ,SIGNXY ,SIGNYZ ,SIGNZX ,THKLY ,
35 7 THK ,SIGY ,ET ,TEMPEL ,DPLANL ,TEMP ,
36 8 SEQ ,INLOC )
37 !=======================================================================
38 ! Implicit types
39 !=======================================================================
40#include "implicit_f.inc"
41 !=======================================================================
42 ! Common
43 !=======================================================================
44 !=======================================================================
45 ! Dummy arguments
46 !=======================================================================
47 INTEGER NEL,NUPARAM,NUVAR,JTHE,INLOC
48 INTEGER ,DIMENSION(NEL), INTENT(IN) :: NGL
49 my_real
50 . TIME,TIMESTEP
51 my_real,DIMENSION(NUPARAM), INTENT(IN) ::
52 . UPARAM
53 my_real,DIMENSION(NEL), INTENT(IN) ::
54 . RHO,TEMPEL,DPLANL,
55 . depsxx,depsyy,depsxy,depsyz,depszx,
56 . sigoxx,sigoyy,sigoxy,sigoyz,sigozx,
57 . gs , thkly
58c
59 my_real ,DIMENSION(NEL), INTENT(OUT) ::
60 . soundsp,sigy,et,
61 . signxx,signyy,signxy,signyz,signzx
62c
63 my_real ,DIMENSION(NEL), INTENT(INOUT) ::
64 . pla,dpla,epsd,off,thk,temp,seq
65 my_real ,DIMENSION(NEL,NUVAR), INTENT(INOUT) ::
66 . uvar
67c
68 !=======================================================================
69 ! Local Variables
70 !=======================================================================
71 INTEGER I,II,NINDX,INDEX(NEL),NICE
72c
73 my_real
74 . young,bulk,lam,g,g2,nu,cdr,kdr,hard,yld0,qvoce,bvoce,jcc,
75 . epsp0,mtemp,tini,tref,eta,cp,dpis,dpad,asrate,afiltr,a11,a12,
76 . nnu, normsig
77 my_real
78 . sigdr2,yld2i,
79 . dtemp,dpdt,dlam,ddep
80 my_real
81 . dpxx,dpyy,dpxy,dpzz,dsdrdj2,dsdrdj3,
82 . dj3dsxx,dj3dsyy,dj3dsxy,dj3dszz,
83 .
84 .
85 . normxx,normyy,normxy,normzz,
86 . denom,dphi,sig_dfdsig,dfdsig2,
87 . dphi_dsig,dphi_dyld,dphi_dfdr,dpdt_nl,
88 . dyld_dpla,dyld_dtemp,dtemp_dlam,dlam_nl
89c
90 my_real, DIMENSION(NEL) ::
91 . dsigxx,dsigyy,dsigxy,trsig,
92 . sxx,syy,sxy,szz,sigm,j2,j3,sigdr,yld,weitemp,
93 . hardp,fhard,frate,ftherm,fdr,phi0,phi,dpla_dlam,
94 . dezz
95 !=======================================================================
96 ! DRUCKER - VOCE - JOHNSON-COOK MATERIAL LAW
97 !=======================================================================
98 !UVAR(1) PHI YIELD FUNCTION VALUE
99 !UVAR(2) EPSD PLASTIC STRAIN RATE SAVED FOR FILTERING
100 !=======================================================================
101c
102 !=======================================================================
103 ! - INITIALISATION OF COMPUTATION ON TIME STEP
104 !=======================================================================
105 ! Recovering model parameters
106 ! Elastic parameters
107 young = uparam(1) ! Young modulus
108 bulk = uparam(2) ! Bulk modulus
109 g = uparam(3) ! Shear modulus
110 g2 = uparam(4) ! 2*Shear modulus
111 lam = uparam(5) ! Lambda Hooke parameter
112 nu = uparam(6) ! Poisson ration
113 nnu = uparam(7) ! NU/(1-NU)
114 a11 = uparam(9) ! Diagonal term, elastic matrix in plane stress
115 a12 = uparam(10) ! Non-diagonal term, elastic matrix in plane stress
116 ! Plastic criterion and hardening parameters [Drucker, 1948]
117 nice = nint(uparam(11)) ! Choice of the Nice method
118 cdr = uparam(12) ! Drucker coefficient
119 kdr = uparam(13) ! Drucker 1/K coefficient
120 tini = uparam(14) ! Initial temperature
121 hard = uparam(15) ! Linear hardening
122 yld0 = uparam(16) ! Initial yield stress
123 qvoce = uparam(17) ! 1st Voce parameter
124 bvoce = uparam(18) ! 2nd Voce parameter
125 ! Strain-rate dependence parameters
126 jcc = uparam(20) ! Johnson-Cook strain rate coefficient
127 epsp0 = uparam(21) ! Johnson-Cook reference strain rate
128 ! Thermal softening and self-heating parameters
129 mtemp = uparam(22) ! Thermal softening slope
130 tref = uparam(23) ! Reference temperature
131 eta = uparam(24) ! Taylor-Quinney coefficient
132 cp = uparam(25) ! Thermal mass capacity
133 dpis = uparam(26) ! Isothermal plastic strain rate
134 dpad = uparam(27) ! Adiabatic plastic strain rate
135 ! Plastic strain-rate filtering parameters
136 asrate = uparam(28) ! Plastic strain rate filtering frequency
137 afiltr = asrate*timestep/(asrate*timestep + one)
138c
139 ! Recovering internal variables
140 DO i=1,nel
141 IF (off(i) < em01) off(i) = zero
142 IF (off(i) < one) off(i) = off(i)*four_over_5
143 ! User inputs
144 phi0(i) = uvar(i,1) ! Previous yield function value
145 ! Standard inputs
146 dpla(i) = zero ! Initialization of the plastic strain increment
147 et(i) = one ! Initialization of hourglass coefficient
148 hardp(i) = zero ! Initialization of hardening modulus
149 ENDDO
150!
151 epsd(1:nel) = uvar(1:nel,2)
152!
153 ! Initialization of temperature and self-heating weight factor
154 IF (time == zero) THEN
155 temp(1:nel) = tini
156 ENDIF
157 IF (cp > zero) THEN
158 IF (jthe == 0) THEN
159 IF (inloc == 0) THEN
160 DO i=1,nel
161 ! Computation of the weight factor
162 IF (epsd(i) < dpis) THEN
163 weitemp(i) = zero
164 ELSEIF (epsd(i) > dpad) THEN
165 weitemp(i) = one
166 ELSE
167 weitemp(i) = ((epsd(i)-dpis)**2 )
168 . * (three*dpad - two*epsd(i) - dpis)
169 . / ((dpad-dpis)**3)
170 ENDIF
171 ENDDO
172 ENDIF
173 ELSE
174 temp(1:nel) = tempel(1:nel)
175 ENDIF
176 ENDIF
177c
178 ! Computation of the initial yield stress
179 DO i = 1,nel
180 ! a) - hardening law
181 fhard(i) = yld0 + hard*pla(i) + qvoce*(one-exp(-bvoce*pla(i)))
182 ! b) - Correction factor for strain-rate dependence (Johnson-Cook)
183 frate(i) = one
184 IF (epsd(i) > epsp0) frate(i) = frate(i) + jcc*log(epsd(i)/epsp0)
185 ! c) - Correction factor for thermal softening
186 ftherm(i) = one
187 IF (cp > zero) ftherm(i) = one - mtemp * (temp(i) - tref)
188 ! d) - Computation of the yield function and value check
189 yld(i) = fhard(i)*frate(i)*ftherm(i)
190 ! e) - Checking values
191 yld(i) = max(em10, yld(i))
192 ENDDO
193c
194 !========================================================================
195 ! - COMPUTATION OF TRIAL VALUES
196 !========================================================================
197 DO i=1,nel
198c
199 ! Computation of the trial stress tensor
200 signxx(i) = sigoxx(i) + a11*depsxx(i) + a12*depsyy(i)
201 signyy(i) = sigoyy(i) + a11*depsyy(i) + a12*depsxx(i)
202 signxy(i) = sigoxy(i) + depsxy(i)*g
203 signyz(i) = sigoyz(i) + depsyz(i)*gs(i)
204 signzx(i) = sigozx(i) + depszx(i)*gs(i)
205 ! computation of the trace of the trial stress tensor
206 trsig(i) = signxx(i) + signyy(i)
207 sigm(i) = -trsig(i) * third
208 ! Computation of the deviatoric trial stress tensor
209 sxx(i) = signxx(i) + sigm(i)
210 syy(i) = signyy(i) + sigm(i)
211 szz(i) = sigm(i)
212 sxy(i) = signxy(i)
213 dezz(i) = -nnu*depsxx(i) - nnu*depsyy(i)
214 ! Second deviatoric invariant
215 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 ) + sxy(i)**2
216 ! Third deviatoric invariant
217 j3(i) = sxx(i)*syy(i)*szz(i) - szz(i)*sxy(i)**2
218 ! Drucker equivalent stress
219 fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
220 ! Checking equivalent stress values
221 IF (fdr(i) > zero) THEN
222 sigdr(i) = kdr * exp((one/six)*log(fdr(i))) ! FDR(I)**(1/6)
223 ELSE
224 sigdr(i) = zero
225 ENDIF
226 ENDDO
227c
228 !========================================================================
229 ! - COMPUTATION OF YIELD FONCTION
230 !========================================================================
231 phi(1:nel) = (sigdr(1:nel) / yld(1:nel))**2 - one
232c
233 ! Checking plastic behavior for all elements
234 nindx = 0
235 DO i=1,nel
236 IF (phi(i) >= zero .AND. off(i) == one) THEN
237 nindx=nindx+1
238 index(nindx)=i
239 ENDIF
240 ENDDO
241c
242 !====================================================================
243 ! - PLASTIC CORRECTION WITH NICE - EXPLICIT METHOD
244 !====================================================================
245#include "vectorize.inc"
246 DO ii=1,nindx
247c
248 ! Number of the element with plastic behaviour
249 i = index(ii)
250c
251 ! Computation of the trial stress increment
252 dsigxx(i) = a11*depsxx(i) + a12*depsyy(i)
253 dsigyy(i) = a11*depsyy(i) + a12*depsxx(i)
254 dsigxy(i) = depsxy(i)*g
255 dlam = zero
256c
257 ! Computation of Drucker equivalent stress of the previous stress tensor
258 trsig(i) = sigoxx(i) + sigoyy(i)
259 sigm(i) = -trsig(i) * third
260 sxx(i) = sigoxx(i) + sigm(i)
261 syy(i) = sigoyy(i) + sigm(i)
262 szz(i) = sigm(i)
263 sxy(i) = sigoxy(i)
264 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2) + sxy(i)**2
265 j3(i) = sxx(i)*syy(i)*szz(i) - szz(i)*sxy(i)**2
266 fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
267 sigdr(i) = kdr * exp((one/six)*log(fdr(i)))
268c
269 ! Note : in this part, the purpose is to compute in one iteration
270 ! a plastic multiplier allowing to update internal variables to satisfy
271 ! the consistency condition.
272 ! Its expression is : DLAMBDA = - (PHI_OLD + DPHI) / DPHI_DLAMBDA
273 ! -> PHI_OLD : old value of yield function (known)
274 ! -> DPHI : yield function prediction (to compute)
275 ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
276 ! into account of internal variables kinetic :
277 ! plasticity, temperature and damage (to compute)
278c
279 ! 1 - Computation of DPHI_DSIG the normal to the yield surface
280 !-------------------------------------------------------------
281c
282 ! Derivative with respect to the equivalent stress
283 yld2i = one/(yld(i)**2)
284 dphi_dsig = two*sigdr(i)*yld2i
285c
286 ! Computation of the Eulerian norm of the stress tensor
287 normsig = sqrt(sigoxx(i)*sigoxx(i)
288 . + sigoyy(i)*sigoyy(i)
289 . + two*sigoxy(i)*sigoxy(i))
290 normsig = max(normsig,one)
291c
292 ! Derivative with respect to Fdr
293 fdr(i) = (j2(i)/(normsig**2))**3 - cdr*((j3(i)/(normsig**3))**2)
294 dphi_dfdr = dphi_dsig*kdr*(one/six)*exp(-(five/six)*log(fdr(i)))
295 dsdrdj2 = dphi_dfdr*three*(j2(i)/(normsig**2))**2
296 dsdrdj3 = -dphi_dfdr*two*cdr*(j3(i)/(normsig**3))
297 ! dJ3/dSig
298 dj3dsxx = two_third*(syy(i)*szz(i))/(normsig**2)
299 . - third*(sxx(i)*szz(i))/(normsig**2)
300 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
301 dj3dsyy = - third*(syy(i)*szz(i))/(normsig**2)
302 . + two_third*(sxx(i)*szz(i))/(normsig**2)
303 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
304 dj3dszz = - third*(syy(i)*szz(i))/(normsig**2)
305 . - third*(sxx(i)*szz(i))/(normsig**2)
306 . + two_third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
307 dj3dsxy = two*(sxx(i)*sxy(i) + sxy(i)*syy(i))/(normsig**2)
308 ! dPhi/dSig
309 normxx = dsdrdj2*sxx(i)/normsig + dsdrdj3*dj3dsxx
310 normyy = dsdrdj2*syy(i)/normsig + dsdrdj3*dj3dsyy
311 normzz = dsdrdj2*szz(i)/normsig + dsdrdj3*dj3dszz
312 normxy = two*dsdrdj2*sxy(i)/normsig + dsdrdj3*dj3dsxy
313c
314 ! Restoring previous value of the yield function
315 phi(i) = phi0(i)
316c
317 ! Computation of yield surface trial increment DPHI
318 dphi = normxx * dsigxx(i)
319 . + normyy * dsigyy(i)
320 . + normxy * dsigxy(i)
321c
322 ! 2 - Computation of DPHI_DLAMBDA
323 !---------------------------------------------------------
324c
325 ! a) Derivative with respect stress increments tensor DSIG
326 ! --------------------------------------------------------
327 dfdsig2 = normxx * (a11*normxx + a12*normyy)
328 . + normyy * (a11*normyy + a12*normxx)
329 . + normxy * normxy * g
330c
331 ! b) derivatives with respect to plastic strain p
332 ! ------------------------------------------------
333c
334 ! i) Derivative of the yield stress with respect to plastic strain dYLD / dPLA
335 ! ----------------------------------------------------------------------------
336 hardp(i) = hard + qvoce*bvoce*exp(-bvoce*pla(i))
337 dyld_dpla = hardp(i)*frate(i)*ftherm(i)
338c
339 ! ii) Derivative of dPLA with respect to DLAM
340 ! -------------------------------------------
341 sig_dfdsig = sigoxx(i) * normxx
342 . + sigoyy(i) * normyy
343 . + sigoxy(i) * normxy
344 dpla_dlam(i) = sig_dfdsig / yld(i)
345c
346 ! c) Derivatives with respect to the temperature TEMP
347 ! ---------------------------------------------------
348 IF (jthe == 0 .AND. cp > zero .AND. inloc == 0) THEN
349 ! i) Derivative of the yield stress with respect to temperature dYLD / dTEMP
350 ! ---------------------------------------------------------------------------
351 dyld_dtemp = -fhard(i)*frate(i)*mtemp
352 ! ii) Derivative of the temperature TEMP with respect to DLAM
353 ! -----------------------------------------------------------
354 dtemp_dlam = weitemp(i)*(eta/(rho(i)*cp))*sig_dfdsig
355 ELSE
356 dyld_dtemp = zero
357 dtemp_dlam = zero
358 ENDIF
359c
360 ! d) Derivative with respect to the yield stress
361 ! ----------------------------------------------
362 sigdr2 = sigdr(i)**2
363 dphi_dyld = -two*sigdr2/(yld(i)**3)
364c
365 ! e) Derivative of PHI with respect to DLAM ( = -DENOM)
366 denom = dfdsig2 - (dphi_dyld * dyld_dpla * dpla_dlam(i))
367 IF (jthe == 0 .AND. cp > zero .AND. inloc == 0) THEN
368 denom = denom - (dphi_dyld * dyld_dtemp * dtemp_dlam)
369 ENDIF
370 denom = sign(max(abs(denom),em20) ,denom)
371c
372 ! 3 - Computation of plastic multiplier and variables update
373 !----------------------------------------------------------
374c
375 ! Computation of the plastic multiplier
376 dlam = (dphi + phi(i)) / denom
377 dlam = max(dlam, zero)
378c
379 ! Plastic strains tensor update
380 dpxx = dlam * normxx
381 dpyy = dlam * normyy
382 dpzz = dlam * normzz
383 dpxy = dlam * normxy
384c
385 ! Elasto-plastic stresses update
386 signxx(i) = sigoxx(i) + dsigxx(i) - (a11*dpxx + a12*dpyy)
387 signyy(i) = sigoyy(i) + dsigyy(i) - (a11*dpyy + a12*dpxx)
388 signxy(i) = sigoxy(i) + dsigxy(i) - dpxy*g
389 trsig(i) = signxx(i) + signyy(i)
390 sigm(i) = -trsig(i) * third
391 sxx(i) = signxx(i) + sigm(i)
392 syy(i) = signyy(i) + sigm(i)
393 szz(i) = sigm(i)
394 sxy(i) = signxy(i)
395c
396 ! Cumulated plastic strain and strain rate update
397 ddep = (dlam/yld(i))*sig_dfdsig
398 dpla(i) = dpla(i) + ddep
399 pla(i) = pla(i) + ddep
400c
401 ! Drucker equivalent stress update
402 j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 ) + sxy(i)**2
403 j3(i) = sxx(i) * syy(i) * szz(i) - szz(i) * sxy(i)**2
404 fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
405 sigdr(i) = kdr * exp((one/six)*log(fdr(i)))
406c
407 ! Temperature update
408 IF (jthe == 0 .AND. cp > zero .AND. inloc == 0) THEN
409 dtemp = weitemp(i)*yld(i)*ddep*eta/(rho(i)*cp)
410 temp(i) = temp(i) + dtemp
411 ftherm(i) = one - mtemp* (temp(i) - tref)
412 ENDIF
413c
414 ! Hardening law update
415 fhard(i) = yld0 + hard * pla(i) + qvoce*(one-exp(-bvoce*pla(i)))
416c
417 ! Yield stress update
418 yld(i) = fhard(i) * frate(i) * ftherm(i)
419 yld(i) = max(yld(i), em10)
420c
421 ! Yield function value update
422 sigdr2 = sigdr(i)**2
423 yld2i = one / yld(i)**2
424 phi(i) = sigdr2 * yld2i - one
425c
426 ! Transverse strain update
427 IF (inloc == 0) THEN
428 dezz(i) = dezz(i) + nnu*dpxx + nnu*dpyy + dpzz
429 ENDIF
430 ENDDO
431 !===================================================================
432 ! - END OF PLASTIC CORRECTION WITH NICE - EXPLICIT METHOD
433 !===================================================================
434c
435 ! Storing new values
436 DO i=1,nel
437 ! USR Outputs & coefficient for hourglass
438 IF (dpla(i) > zero) THEN
439 uvar(i,1) = phi(i) ! Yield function value
440 et(i) = hardp(i)*frate(i) / (hardp(i)*frate(i) + young)
441 ELSE
442 uvar(i,1) = zero
443 et(i) = one
444 ENDIF
445 seq(i) = sigdr(i)
446 ! Plastic strain-rate (filtered)
447 dpdt = dpla(i) / max(em20,timestep)
448 epsd(i) = afiltr * dpdt + (one - afiltr) * epsd(i)
449 uvar(i,2) = epsd(i)
450 ! Computation of the sound speed
451 soundsp(i) = sqrt(a11/rho(i))
452 ! Storing the yield stress
453 sigy(i) = yld(i)
454 ! Non-local thickness variation + temperature
455 IF ((inloc > 0).AND.(off(i) == one).AND.(dplanl(i)>zero)) THEN
456 yld2i = one/(yld(i)**2)
457 dphi_dsig = two*sigdr(i)*yld2i
458 normsig = sqrt(signxx(i)*signxx(i)
459 . + signyy(i)*signyy(i)
460 . + two*signxy(i)*signxy(i))
461 normsig = max(normsig,one)
462 fdr(i) = (j2(i)/(normsig**2))**3 - cdr*((j3(i)/(normsig**3))**2)
463 IF (fdr(i) > zero) THEN
464 dphi_dfdr = dphi_dsig*kdr*(one/six)*exp(-(five/six)*log(fdr(i)))
465 ELSE
466 dphi_dfdr = zero
467 ENDIF
468 dsdrdj2 = dphi_dfdr*three*(j2(i)/(normsig**2))**2
469 dsdrdj3 = -dphi_dfdr*two*cdr*(j3(i)/(normsig**3))
470 dj3dsxx = two_third*(syy(i)*szz(i))/(normsig**2)
471 . - third*(sxx(i)*szz(i))/(normsig**2)
472 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
473 dj3dsyy = - third*(syy(i)*szz(i))/(normsig**2)
474 . + two_third*(sxx(i)*szz(i))/(normsig**2)
475 . - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
476 dj3dszz = - third*(syy(i)*szz(i))/(normsig**2)
477 . - third*(sxx(i)*szz(i))/(normsig**2)
478 . + two_third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)
479 dj3dsxy = two*(sxx(i)*sxy(i) + sxy(i)*syy(i))/(normsig**2)
480 normxx = dsdrdj2*sxx(i)/normsig + dsdrdj3*dj3dsxx
481 normyy = dsdrdj2*syy(i)/normsig + dsdrdj3*dj3dsyy
482 normzz = dsdrdj2*szz(i)/normsig + dsdrdj3*dj3dszz
483 normxy = two*dsdrdj2*sxy(i)/normsig + dsdrdj3*dj3dsxy
484 sig_dfdsig = signxx(i) * normxx
485 . + signyy(i) * normyy
486 . + signxy(i) * normxy
487 IF (sig_dfdsig > em01) THEN
488 dlam_nl = (yld(i)*dplanl(i))/sig_dfdsig
489 dezz(i) = dezz(i) + nnu*(dlam_nl*normxx)
490 . + nnu*(dlam_nl*normyy)
491 . + dlam_nl*normzz
492 ENDIF
493 IF (cp > zero .AND. jthe == 0) THEN
494 ! Computation of the weight factor
495 dpdt_nl = dplanl(i)/max(timestep,em20)
496 IF (dpdt_nl < dpis) THEN
497 weitemp(i) = zero
498 ELSEIF (dpdt_nl > dpad) THEN
499 weitemp(i) = one
500 ELSE
501 weitemp(i) = ((dpdt_nl-dpis)**2 )
502 . * (three*dpad - two*dpdt_nl - dpis)
503 . / ((dpad-dpis)**3)
504 ENDIF
505 dtemp = weitemp(i)*dplanl(i)*yld(i)*eta/(rho(i)*cp)
506 temp(i) = temp(i) + dtemp
507 ENDIF
508 ENDIF
509 ! Computation of the thickness variation
510 thk(i) = thk(i) + dezz(i)*thkly(i)*off(i)
511 ENDDO
512c
513 END
the
end diagonal values have been computed in the(sparse) matrix id.SOL
max
#define max(a, b)
Definition macros.h:21
mat104c_nodam_nice
subroutine mat104c_nodam_nice(nel, ngl, nuparam, nuvar, time, timestep, uparam, uvar, jthe, off, gs, rho, pla, dpla, epsd, soundsp, depsxx, depsyy, depsxy, depsyz, depszx, sigoxx, sigoyy, sigoxy, sigoyz, sigozx, signxx, signyy, signxy, signyz, signzx, thkly, thk, sigy, et, tempel, dplanl, temp, seq, inloc)
Definition mat104c_nodam_nice.F:37