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