OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps121.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!|| sigeps121 ../engine/source/materials/mat/mat121/sigeps121.F
25!||--- called by ------------------------------------------------------
26!|| mulaw ../engine/source/materials/mat_share/mulaw.F90
27!||--- calls -----------------------------------------------------
28!|| mat121_newton ../engine/source/materials/mat/mat121/mat121_newton.F
29!|| mat121_nice ../engine/source/materials/mat/mat121/mat121_nice.F
30!|| vinter2 ../engine/source/tools/curve/vinter.F
31!||--- uses -----------------------------------------------------
32!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
33!||====================================================================
34 SUBROUTINE sigeps121(
35 1 NEL ,NGL ,NUPARAM ,NUVAR ,NFUNC ,IFUNC ,NPF ,
36 2 TF ,TIMESTEP,TIME ,UPARAM ,UVAR ,RHO ,PLA ,
37 3 DPLA ,SOUNDSP ,EPSD ,OFF ,LOFF ,
38 4 DEPSXX ,DEPSYY ,DEPSZZ ,DEPSXY ,DEPSYZ ,DEPSZX ,
39 5 EPSPXX ,EPSPYY ,EPSPZZ ,EPSPXY ,EPSPYZ ,EPSPZX ,
40 6 SIGOXX ,SIGOYY ,SIGOZZ ,SIGOXY ,SIGOYZ ,SIGOZX ,
41 7 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
42 8 SIGY ,ET ,VARNL ,INLOC ,DT ,
43 9 IPG ,NPG ,ELBUF_TAB)
44C-----------------------------------------------
45C M o d u l e s
46C-----------------------------------------------
47 USE elbufdef_mod
48C-----------------------------------------------
49C I m p l i c i t T y p e s
50C-----------------------------------------------
51#include "implicit_f.inc"
52C-----------------------------------------------
53C C O M M O N
54C-----------------------------------------------
55#include "scr17_c.inc"
56#include "com08_c.inc"
57#include "units_c.inc"
58#include "comlock.inc"
59C-----------------------------------------------
60C D u m m y A r g u m e n t s
61C-----------------------------------------------
62 INTEGER NEL,NUPARAM,NUVAR,NFUNC,INLOC,
63 . IPG,NPG,NPF(*),NGL(NEL),
64 . IFUNC(NFUNC)
65 my_real
66 . TIME,TIMESTEP,TF(*),UPARAM(NUPARAM)
67 my_real,DIMENSION(NEL), INTENT(IN) ::
68 . RHO,DT,
69 . DEPSXX,DEPSYY,DEPSZZ,DEPSXY,DEPSYZ,DEPSZX,
70 . EPSPXX,EPSPYY,EPSPZZ,EPSPXY,EPSPYZ,EPSPZX,
71 . sigoxx,sigoyy,sigozz,sigoxy,sigoyz,sigozx
72 my_real ,DIMENSION(NEL), INTENT(OUT) ::
73 . soundsp,signxx,signyy,signzz,signxy,signyz,signzx
74 my_real,DIMENSION(NEL) ::
75 . sigy,et
76 my_real ,DIMENSION(NEL), INTENT(INOUT) ::
77 . pla,dpla,epsd,varnl,loff,off
78 my_real ,DIMENSION(NEL,NUVAR), INTENT(INOUT) ::
79 . uvar
80 TYPE (ELBUF_STRUCT_), TARGET :: ELBUF_TAB
81C-----------------------------------------------
82C L o c a l V a r i a b l e s
83C-----------------------------------------------
84 INTEGER I,J,IRES,Ifail,NINDX,NINDX2,INDX(NEL),INDX2(NEL),
85 . ipos(nel),iad(nel),ilen(nel),ir,is,it
86 my_real dtmin,xscale_fail,yscale_fail,s1,s2,q,s11,
87 . s22,r_inter,dfdepsd(nel),fail(nel),seq(nel),
88 . i1,i2,i3,r,psi,s33
89C=======================================================================
90c
91 ires = nint(uparam(11)) ! Plastic projection method
92 ! = 1 => Nice method
93 ! = 2 => Newton-iteration method
94 ifail = nint(uparam(13)) ! Failure criterion flag
95 ! = 0 => Von Mises stress
96 ! = 1 => Plastic strain
97 ! = 2 => Maximum princ. stress +
98 ! abs(Minimum princ. stress)
99 ! = 3 => Maximum princ. stress
100 dtmin = uparam(15) ! Minimal timestep for element deletion
101 xscale_fail = uparam(22) ! Strain-rate scale factor for failure criterion function
102 yscale_fail = uparam(23) ! Ordinate scale factor for failure criterion function
103c--------------------------
104 SELECT CASE (ires)
105c
106 CASE(1)
107c
108 CALL mat121_nice(
109 1 nel ,ngl ,nuparam ,nuvar ,nfunc ,ifunc ,npf ,
110 2 tf ,timestep,time ,uparam ,uvar ,rho ,pla ,
111 3 dpla ,soundsp ,epsd ,off ,
112 4 depsxx ,depsyy ,depszz ,depsxy ,depsyz ,depszx ,
113 5 epspxx ,epspyy ,epspzz ,epspxy ,epspyz ,epspzx ,
114 6 sigoxx ,sigoyy ,sigozz ,sigoxy ,sigoyz ,sigozx ,
115 7 signxx ,signyy ,signzz ,signxy ,signyz ,signzx ,
116 8 sigy ,et ,seq )
117c
118 CASE(2)
119c
120 CALL mat121_newton(
121 1 nel ,ngl ,nuparam ,nuvar ,nfunc ,ifunc ,npf ,
122 2 tf ,timestep,time ,uparam ,uvar ,rho ,pla ,
123 3 dpla ,soundsp ,epsd ,off ,
124 4 depsxx ,depsyy ,depszz ,depsxy ,depsyz ,depszx ,
125 5 epspxx ,epspyy ,epspzz ,epspxy ,epspyz ,epspzx ,
126 6 sigoxx ,sigoyy ,sigozz ,sigoxy ,sigoyz ,sigozx ,
127 7 signxx ,signyy ,signzz ,signxy ,signyz ,signzx ,
128 8 sigy ,et ,seq )
129c
130 END SELECT
131c
132 !--------------------------------------------------------------------
133 ! Failure
134 !--------------------------------------------------------------------
135c
136 ! Compute failure criterion value
137 IF (ifunc(4) > 0) THEN
138 ipos(1:nel) = 1
139 iad(1:nel) = npf(ifunc(4)) / 2 + 1
140 ilen(1:nel) = npf(ifunc(4)+1) / 2 - iad(1:nel) - ipos(1:nel)
141 CALL vinter2(tf,iad,ipos,ilen,nel,epsd/xscale_fail,dfdepsd,fail)
142 fail(1:nel) = yscale_fail*fail(1:nel)
143 ENDIF
144c
145 ! Checking elements deletion criteria
146 nindx = 0
147 nindx2 = 0
148 IF (dtmin > zero .OR. ifunc(4) > 0) THEN
149 DO i = 1,nel
150c
151 ! Minimum timestep
152 IF ((dt(i) > em20).AND.(dt(i) < dtmin).AND.(off(i) == one)) THEN
153 off(i) = zero
154 nindx = nindx+1
155 indx(nindx) = i
156 ENDIF
157c
158 ! Failure criterion
159 IF ((ifunc(4) > 0).AND.(off(i) == one)) THEN
160 ! Failure indicator evolution
161 IF (loff(i) < em01) loff(i) = zero
162 IF (loff(i) < one) loff(i) = loff(i)*four_over_5
163c
164 ! Fully integrated solid elements
165 IF (npg > 1) THEN
166 ! Checking full failure of the element
167 IF (ipg == npg) THEN
168 ! Initialization of OFFG
169 off(i) = zero
170 ! Loop over integration points
171 DO ir = 1, elbuf_tab%NPTR
172 DO is = 1, elbuf_tab%NPTS
173 DO it = 1, elbuf_tab%NPTT
174 !If one integration points is not fully broken, the brick remains
175 IF (elbuf_tab%BUFLY(1)%LBUF(ir,is,it)%OFF(i)>zero) off(i) = one
176 ENDDO
177 ENDDO
178 ENDDO
179 ENDIF
180 ! Under-integrated solid element
181 ELSE
182 !initialization for checking complete failure of the shell(all integration points)
183 IF (ipg == 1) THEN
184 off(i) = zero
185 ENDIF
186 !If one integration points is not fully broken, the brick remains
187 IF (loff(i)>zero) off(i) = one
188 ENDIF
189 ! Check integration point failure
190 IF (loff(i) == one) THEN
191 ! -> Von Mises stress
192 IF (ifail == 0) THEN
193 IF (seq(i) >= fail(i)) loff(i) = four_over_5
194 ! -> Plastic strain
195 ELSEIF (ifail == 1) THEN
196 IF (pla(i) >= fail(i)) loff(i) = four_over_5
197 ! -> Maximum principal stress and absolute value of minimum principal stress
198 ELSEIF (ifail == 2) THEN
199 ! Computing the principal stresses
200 i1 = signxx(i)+signyy(i)+signzz(i)
201 i2 = signxx(i)*signyy(i)+signyy(i)*signzz(i)+signzz(i)*signxx(i)-
202 . signxy(i)*signxy(i)-signzx(i)*signzx(i)-signyz(i)*signyz(i)
203 i3 = signxx(i)*signyy(i)*signzz(i)-signxx(i)*signyz(i)*signyz(i)-
204 . signyy(i)*signzx(i)*signzx(i)-signzz(i)*signxy(i)*signxy(i)+
205 . two*signxy(i)*signzx(i)*signyz(i)
206 q = (three*i2 - i1*i1)/nine
207 r = (two*i1*i1*i1-nine*i1*i2+twenty7*i3)/cinquante4 ! (2*I3^3-9*I1*I2+27*I3)/54
208 r_inter = min(r/sqrt(max(em20,(-q**3))),one)
209 psi = acos(max(r_inter,-one))
210 s11 = two*sqrt(-q)*cos(psi/three)+third*i1
211 s22 = two*sqrt(-q)*cos((psi+two*pi)/three)+third*i1
212 s33 = two*sqrt(-q)*cos((psi+four*pi)/three)+third*i1
213 ! Sorting principal strains
214 IF (s11 < s22) THEN
215 r_inter = s11
216 s11 = s22
217 s22 = r_inter
218 ENDIF
219 IF (s22 < s33)THEN
220 r_inter = s22
221 s22 = s33
222 s33 = r_inter
223 ENDIF
224 IF (s11 < s22)THEN
225 r_inter = s11
226 s11 = s22
227 s22 = r_inter
228 ENDIF
229 IF ((s11>=fail(i)).OR.(abs(s33)>=fail(i))) loff(i) = four_over_5
230 ! -> Maximum principal stress
231 ELSEIF (ifail == 3) THEN
232 ! Computing the principal stresses
233 i1 = signxx(i)+signyy(i)+signzz(i)
234 i2 = signxx(i)*signyy(i)+signyy(i)*signzz(i)+signzz(i)*signxx(i)-
235 . signxy(i)*signxy(i)-signzx(i)*signzx(i)-signyz(i)*signyz(i)
236 i3 = signxx(i)*signyy(i)*signzz(i)-signxx(i)*signyz(i)*signyz(i)-
237 . signyy(i)*signzx(i)*signzx(i)-signzz(i)*signxy(i)*signxy(i)+
238 . two*signxy(i)*signzx(i)*signyz(i)
239 q = (three*i2 - i1*i1)/nine
240 r = (two*i1*i1*i1-nine*i1*i2+twenty7*i3)/cinquante4 ! (2*I3^3-9*I1*I2+27*I3)/54
241 r_inter = min(r/sqrt(max(em20,(-q**3))),one)
242 psi = acos(max(r_inter,-one))
243 s11 = two*sqrt(-q)*cos(psi/three)+third*i1
244 s22 = two*sqrt(-q)*cos((psi+two*pi)/three)+third*i1
245 s33 = two*sqrt(-q)*cos((psi+four*pi)/three)+third*i1
246 ! Sorting principal strains
247 IF (s11 < s22) THEN
248 r_inter = s11
249 s11 = s22
250 s22 = r_inter
251 ENDIF
252 IF (s22 < s33)THEN
253 r_inter = s22
254 s22 = s33
255 s33 = r_inter
256 ENDIF
257 IF (s11 < s22)THEN
258 r_inter = s11
259 s11 = s22
260 s22 = r_inter
261 ENDIF
262 IF (s11>=fail(i)) loff(i) = four_over_5
263 ENDIF
264 !Integration point failure
265 IF (loff(i) == four_over_5) THEN
266 nindx2 = nindx2+1
267 indx2(nindx2) = i
268 ENDIF
269 ENDIF
270 ENDIF
271 ENDDO
272 ENDIF
273c
274 ! Printout timestep element deletion
275 IF (nindx>0) THEN
276 DO j=1,nindx
277#include "lockon.inc"
278 WRITE(iout, 1000) ngl(indx(j))
279 WRITE(istdo,1100) ngl(indx(j)),tt
280#include "lockoff.inc"
281 ENDDO
282 ENDIF
283c
284 ! Printout failure
285 IF (nindx2>0) THEN
286 DO j=1,nindx2
287#include "lockon.inc"
288 WRITE(iout, 2000) ngl(indx2(j)),ipg
289 WRITE(istdo,2100) ngl(indx2(j)),ipg,tt
290#include "lockoff.inc"
291 ENDDO
292 ENDIF
293c
294 1000 FORMAT(1x,'MINIMUM TIMESTEP (PLAS_RATE) REACHED, DELETED SOLID ELEMENT ',i10)
295 1100 FORMAT(1x,'MINIMUM TIMESTEP (PLAS_RATE) REACHED, DELETED SOLID ELEMENT ',i10,1x,'AT TIME :',1pe12.4)
296 2000 FORMAT(1x,'FAILURE (PLAS_RATE) IN SOLID ELEMENT ',i10,1x,',GAUSS PT',i2)
297 2100 FORMAT(1x,'FAILURE (PLAS_RATE) IN SOLID ELEMENT ',i10,1x,',GAUSS PT',i2,1x,'AT TIME :',1pe12.4)
298c
299c-----------
300 RETURN
301 END
the
end diagonal values have been computed in the(sparse) matrix id.SOL
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
mat121_newton
subroutine mat121_newton(nel, ngl, nuparam, nuvar, nfunc, ifunc, npf, tf, timestep, time, uparam, uvar, rho, pla, dpla, soundsp, epsd, off, depsxx, depsyy, depszz, depsxy, depsyz, depszx, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, sigy, et, seq)
Definition mat121_newton.F:39
mat121_nice
subroutine mat121_nice(nel, ngl, nuparam, nuvar, nfunc, ifunc, npf, tf, timestep, time, uparam, uvar, rho, pla, dpla, soundsp, epsd, off, depsxx, depsyy, depszz, depsxy, depsyz, depszx, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, sigy, et, seq)
Definition mat121_nice.F:39
for
for(i8=*sizetab-1;i8 >=0;i8--)
Definition mumps_common.c:91
sigeps121
subroutine sigeps121(nel, ngl, nuparam, nuvar, nfunc, ifunc, npf, tf, timestep, time, uparam, uvar, rho, pla, dpla, soundsp, epsd, off, loff, depsxx, depsyy, depszz, depsxy, depsyz, depszx, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, sigy, et, varnl, inloc, dt, ipg, npg, elbuf_tab)
Definition sigeps121.F:44
vinter2
subroutine vinter2(tf, iad, ipos, ilen, nel0, x, dydx, y)
Definition vinter.F:143