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