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