OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps79.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!|| sigeps79 ../engine/source/materials/mat/mat079/sigeps79.F
25!||--- called by ------------------------------------------------------
26!|| mulaw ../engine/source/materials/mat_share/mulaw.F90
27!||====================================================================
28 SUBROUTINE sigeps79(
29 1 NEL ,NUPARAM ,NUVAR ,TIME ,TIMESTEP ,UPARAM ,
30 2 RHO0 ,RHO ,NGL ,SIGY ,DPLA ,DEFP ,
31 3 DEPSXX ,DEPSYY ,DEPSZZ ,DEPSXY ,DEPSYZ ,DEPSZX ,
32 4 SIGOXX ,SIGOYY ,SIGOZZ ,SIGOXY ,SIGOYZ ,SIGOZX ,
33 5 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
34 6 EPSD ,DMG ,SOUNDSP ,UVAR ,OFF ,AMU ,
35 7 ET )
36C
37C-----------------------------------------------
38C I M P L I C I T T Y P E S
39C-----------------------------------------------
40#include "implicit_f.inc"
41C-----------------------------------------------
42C C O M M O N
43C-----------------------------------------------
44#include "param_c.inc"
45#include "comlock.inc"
46#include "units_c.inc"
47C----------------------------------------------------------------
48C I N P U T A R G U M E N T S
49C----------------------------------------------------------------
50 INTEGER,INTENT(IN) ::
51 . NEL,NUPARAM,NUVAR,NGL(NEL)
52 my_real,INTENT(IN) ::
53 . TIME,TIMESTEP,UPARAM(NUPARAM),RHO0(NEL),RHO(NEL),
54 . DEPSXX(NEL),DEPSYY(NEL),DEPSZZ(NEL),
55 . DEPSXY(NEL),DEPSYZ(NEL),DEPSZX(NEL),
56 . SIGOXX(NEL),SIGOYY(NEL),SIGOZZ(NEL),
57 . sigoxy(nel),sigoyz(nel),sigozx(nel),
58 . epsd(nel),amu(nel)
59C----------------------------------------------------------------
60C O U T P U T A R G U M E N T S
61C----------------------------------------------------------------
62 my_real,INTENT(OUT) ::
63 . sigy(nel),dpla(nel),soundsp(nel),
64 . signxx(nel),signyy(nel),signzz(nel),
65 . signxy(nel),signyz(nel),signzx(nel)
66C----------------------------------------------------------------
67C I N P U T O U T P U T A R G U M E N T S
68C----------------------------------------------------------------
69 my_real,INTENT(INOUT) ::
70 . uvar(nel,nuvar),off(nel),dmg(nel),
71 . defp(nel),et(nel)
72C----------------------------------------------------------------
73C L O C A L V A R I A B L E S
74C----------------------------------------------------------------
75 INTEGER I,J,NINDX,INDX(NEL),IDEL
76 my_real
77 . g , g2 , aa , bb , mm ,
78 . nn , cc , eps0 , sigfmax,
79 . tstar , phel , shel , beta ,
80 . d1 , d2 , epsmax, k1 , k2 , k3
81 my_real
82 . mu(nel),mu2(nel),pold(nel),vm(nel),
83 . deltap(nel),pnew(nel),pstar(nel),scale(nel),
84 . sigyi(nel),sigyf(nel),sigyold(nel),dmg_old(nel)
85 my_real
86 . dav, ce, sigstar, epfail, p1, yield, deltau, dpdmu,
87 . pmin, ratio, j2
88C
89 !=======================================================================
90 ! - INITIALISATION OF COMPUTATION ON TIME STEP
91 !=======================================================================
92 ! Recovering model parameters
93 g = uparam(1)
94 g2 = uparam(2)
95 aa = uparam(3)
96 bb = uparam(4)
97 mm = uparam(5)
98 nn = uparam(6)
99 cc = uparam(7)
100 eps0 = uparam(8)
101 sigfmax = uparam(9)
102 tstar = uparam(10)
103 phel = uparam(11)
104 shel = uparam(12)
105 d1 = uparam(13)
106 d2 = uparam(14)
107 k1 = uparam(15)
108 k2 = uparam(16)
109 k3 = uparam(17)
110 beta = uparam(18)
111 idel = nint(uparam(19))
112 epsmax = uparam(20)
113c
114 ! Recovering internal variables
115 DO i=1,nel
116 IF (off(i) < em01) off(i) = zero
117 IF (off(i) < one) off(i) = off(i)*four_over_5
118 deltap(i) = uvar(i,1)
119 sigyold(i) = uvar(i,2)/shel
120 dmg_old(i) = dmg(i)
121 mu(i) = amu(i)
122 mu2(i) = mu(i)*mu(i)
123 ENDDO
124c
125 !========================================================================
126 ! - COMPUTATION OF ELASTIC DEVIATORIC STRESSES AND EQUIVALENT STRESS
127 !========================================================================
128 DO i=1,nel
129 dav = (depsxx(i)+depsyy(i)+depszz(i))*third
130 pold(i) = -(sigoxx(i)+sigoyy(i)+sigozz(i))*third
131 signxx(i) = sigoxx(i)+pold(i)+g2*(depsxx(i)-dav)
132 signyy(i) = sigoyy(i)+pold(i)+g2*(depsyy(i)-dav)
133 signzz(i) = sigozz(i)+pold(i)+g2*(depszz(i)-dav)
134 signxy(i) = sigoxy(i)+g*depsxy(i)
135 signyz(i) = sigoyz(i)+g*depsyz(i)
136 signzx(i) = sigozx(i)+g*depszx(i)
137 j2 = half*(signxx(i)**2+signyy(i)**2+signzz(i)**2)
138 . +signxy(i)**2+signyz(i)**2+signzx(i)**2
139 vm(i) = sqrt(three*j2)
140 ENDDO
141c
142 !========================================================================
143 ! - COMPUTATION OF THE PRESSURE
144 !========================================================================
145 DO i=1,nel
146 pnew(i) = k1*mu(i) + deltap(i)
147 IF (mu(i) > zero) THEN
148 pnew(i) = pnew(i) + k2*mu2(i) + k3*mu2(i)*mu(i)
149 ELSEIF (idel /= 1) THEN
150 pmin=-tstar*phel*(one-dmg(i))
151 pnew(i)=max(pnew(i),pmin)
152 ENDIF
153 pstar(i) = pnew(i)/phel
154 ENDDO
155c
156 !========================================================================
157 ! - COMPUTATION OF THE YIELD STRESS
158 !========================================================================
159 DO i=1,nel
160 IF (nn == zero) THEN
161 sigyi(i) = aa
162 ELSEIF ((pstar(i)+tstar) > zero) THEN
163 sigyi(i) = aa*(pstar(i)+tstar)**nn
164 ELSE
165 sigyi(i) = zero
166 ENDIF
167C
168 IF (mm == zero) THEN
169 sigyf(i) = bb
170 ELSEIF (pstar(i) > zero) THEN
171 sigyf(i) = bb*(pstar(i))**mm
172 ELSE
173 sigyf(i) = zero
174 ENDIF
175C
176 IF (epsd(i)<=eps0) THEN
177 ce = one
178 ELSE
179 ce = one + cc*log(epsd(i)/eps0)
180 ENDIF
181C
182 sigyi(i) = ce*sigyi(i)
183 sigyf(i) = ce*sigyf(i)
184 sigyf(i) = min(sigyf(i),sigfmax)
185 sigy(i) = (one-dmg(i))*sigyi(i)+dmg(i)*sigyf(i)
186 ENDDO
187c
188 !========================================================================
189 ! - RADIAL RETURN
190 !========================================================================
191 DO i=1,nel
192 IF (off(i) == one) THEN
193 sigstar = vm(i)/shel
194 IF (sigstar < sigy(i)) THEN
195 scale(i) = one
196 ELSEIF (vm(i) > zero) THEN
197 scale(i) = sigy(i)/sigstar
198 ELSE
199 scale(i) = zero
200 ENDIF
201 signxx(i) = scale(i)*signxx(i)
202 signyy(i) = scale(i)*signyy(i)
203 signzz(i) = scale(i)*signzz(i)
204 signxy(i) = scale(i)*signxy(i)
205 signyz(i) = scale(i)*signyz(i)
206 signzx(i) = scale(i)*signzx(i)
207 ENDIF
208 ENDDO
209c
210 !========================================================================
211 ! - UPDATE PLASTIC STRAIN AND DAMAGE
212 !========================================================================
213 nindx = 0
214 indx(1:nel) = 0
215 DO i=1,nel
216 IF (off(i) == one) THEN
217c
218 ! Compute plastic strain at failure
219 IF (d2 == zero) THEN
220 epfail = d1
221 ELSEIF ((pstar(i)+tstar) >= zero) THEN
222 epfail = d1*(pstar(i)+tstar)**d2
223 ELSE
224 epfail = zero
225 ENDIF
226c
227 ! Update plastic strain and damage
228 ! -> When damage parameters are defined, progressive failure + plastic strain
229 IF (epfail > zero) THEN
230 dpla(i) = (one - scale(i))*vm(i)/(three*sqrt(three)*g)
231 defp(i) = defp(i) + dpla(i)
232 dmg(i) = dmg(i) + dpla(i)/epfail
233 dmg(i) = min(dmg(i),one)
234 ! -> When damage parameters are not defined, no plastic strain and instant failure
235 ELSEIF (scale(i)<one) THEN
236 dmg(i) = one
237 ENDIF
238c
239 ! Check element deletion
240 IF (idel == 1) THEN
241 IF ((pstar(i)+tstar) < zero) THEN
242 off(i) = four_over_5
243 nindx = nindx + 1
244 indx(nindx) = i
245 ENDIF
246 ELSEIF (idel == 2) THEN
247 IF (defp(i) > epsmax) THEN
248 off(i) = four_over_5
249 nindx = nindx + 1
250 indx(nindx) = i
251 ENDIF
252 ELSEIF (idel == 3) THEN
253 IF (dmg(i) == one) THEN
254 off(i) = four_over_5
255 nindx = nindx + 1
256 indx(nindx) = i
257 ENDIF
258 ENDIF
259 ENDIF
260 ENDDO
261c
262 !========================================================================
263 ! - COMPUTE PRESSURE INCREMENT
264 !========================================================================
265 DO i=1,nel
266 IF ((dmg(i) > dmg_old(i)).AND.(mu(i) > zero).AND.(off(i) == one)) THEN
267 p1 = k1*mu(i)
268 yield = (one-dmg(i))*sigyi(i)+dmg(i)*sigyf(i)
269 deltau = (sigyold(i)*sigyold(i)-yield*yield)/(six*g)
270 IF (deltau > zero) THEN
271 deltau = deltau*shel*shel
272 deltap(i) = -p1+
273 . sqrt((deltap(i)+p1)**2+two*beta*k1*deltau)
274 ENDIF
275 ENDIF
276 ENDDO
277C
278 !========================================================================
279 ! - UPDATE STRESS TENSOR AND SOUND SPEED
280 !========================================================================
281 DO i=1,nel
282 uvar(i,1) = deltap(i)
283 uvar(i,2) = sigy(i)*shel
284 signxx(i) = signxx(i)-pnew(i)
285 signyy(i) = signyy(i)-pnew(i)
286 signzz(i) = signzz(i)-pnew(i)
287 IF (mu(i) > zero) THEN
288 dpdmu = k1+two*k2*mu(i)+three*k3*mu2(i)
289 ELSE
290 dpdmu = k1
291 ENDIF
292 soundsp(i) = sqrt((dpdmu+four_over_3*g)/rho0(i))
293 et(i) = one
294 ENDDO
295C
296 !========================================================================
297 ! - PRINT OUT ELEMENT DELETION DATA
298 !========================================================================
299 IF (nindx>0)THEN
300 DO j=1,nindx
301#include "lockon.inc"
302 WRITE(iout, 1000) ngl(indx(j)),time
303 WRITE(istdo,1000) ngl(indx(j)),time
304#include "lockoff.inc"
305 ENDDO
306 ENDIF
307C
308 1000 FORMAT(1x,'-- RUPTURE (J-HOLMQUIST) OF SOLID ELEMENT :',i10,' AT TIME :',1pe12.4)
309C
310 END
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
sigeps79
subroutine sigeps79(nel, nuparam, nuvar, time, timestep, uparam, rho0, rho, ngl, sigy, dpla, defp, depsxx, depsyy, depszz, depsxy, depsyz, depszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, epsd, dmg, soundsp, uvar, off, amu, et)
Definition sigeps79.F:36