OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps88c.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!|| sigeps88c ../engine/source/materials/mat/mat088/sigeps88c.F
25!||--- called by ------------------------------------------------------
26!|| mulawc ../engine/source/materials/mat_share/mulawc.F90
27!||--- calls -----------------------------------------------------
28!|| finter ../engine/source/tools/curve/finter.F
29!||====================================================================
30 SUBROUTINE sigeps88c(
31 1 NEL , NUPARAM, NUVAR , NFUNC , IFUNC , NPF ,
32 2 NPT0 , ILAYER , NGL , OFF , ISMSTR, GS ,
33 2 TF , TIME , TIMESTEP, UPARAM, RHO0 ,
34 3 AREA , EINT , THKLYL,
35 5 DEPSXX , DEPSYY , DEPSXY, DEPSYZ, DEPSZX,
36 6 EPSXX , EPSYY , EPSXY , EPSYZ , EPSZX ,
37 7 SIGOXX , SIGOYY , SIGOXY, SIGOYZ, SIGOZX,
38 8 SIGNXX , SIGNYY , SIGNXY, SIGNYZ, SIGNZX,
39 A SOUNDSP, VISCMAX, THKN , UVAR )
40C-----------------------------------------------
41C I M P L I C I T T Y P E S
42C-----------------------------------------------
43#include "implicit_f.inc"
44C-----------------------------------------------
45C C O M M O N
46C-----------------------------------------------
47#include "tabsiz_c.inc"
48C----------------------------------------------------------------
49C I N P U T A R G U M E N T S
50C----------------------------------------------------------------
51 INTEGER NEL,NUPARAM,NUVAR,ISMSTR,NPT0,ILAYER
52 INTEGER, DIMENSION(NEL) :: NGL
53 my_real
54 . UPARAM(NUPARAM)
55 my_real , DIMENSION(NEL,2) :: EINT
56 my_real
57 . TIME,TIMESTEP
58 my_real, DIMENSION(NEL) ::
59 . THKN,THKLYL, RHO0,AREA,GS,
60 . depsxx,depsyy,depsxy,depsyz,depszx,
61 . epsxx ,epsyy ,epsxy ,epsyz ,epszx ,
62 . sigoxx,sigoyy,sigoxy,sigoyz,sigozx
63C----------------------------------------------------------------
64C O U T P U T A R G U M E N T S
65C----------------------------------------------------------------
66 my_real, DIMENSION(NEL) ::
67 . signxx ,signyy ,signxy ,signyz ,signzx,soundsp,viscmax
68C----------------------------------------------------------------
69C I N P U T O U T P U T A R G U M E N T S
70C----------------------------------------------------------------
71 my_real
72 . uvar(nel,nuvar), off(nel)
73C----------------------------------------------------------------
74C VARIABLES FOR FUNCTION INTERPOLATION
75C----------------------------------------------------------------
76 INTEGER NPF(SNPC), NFUNC, IFUNC(NFUNC)
77 my_real FINTER,FINTTE,TF(STF),FINT2V
78 EXTERNAL FINTER,FINTTE
79C----------------------------------------------------------------
80C L O C A L V A R I B L E S
81C----------------------------------------------------------------
82 INTEGER IUNL_FOR,ICASE,I,ITER,NLOAD,NN,J,INDX_L(NEL),INDX_UNL(NEL),
83 . NE_L,NE_UNL,II
84 my_real
85 . rbulk, nu,hys, shape,kt,y3,a11,dd2,dd1,dd3,df1,df2,df3,
86 . dflam3, emax, fac,y33,df, p,invr,yfac(nfunc),rate(nfunc),yfac_unl,
87 . scale,g_ini,dam
88 my_real, DIMENSION(NEL) :: t1,t2,t3,epsz,trav,rootv,f1,f2,f3,gmax,
89 . rv,dezz,epszz,ecurent,deint,eps_eq,deps_eq
90 .
91 my_real :: ee(nel,3),evv(nel,3),ev(nel,3), eigv(nel,3,3)
92C=======================================================================
93C material parameters
94 rbulk = uparam(1)
95 nu = uparam(2)
96 gs = uparam(3)
97 nload = int(uparam(4))
98 iunl_for = int(uparam(5))
99 hys = uparam(6)
100 shape = uparam(7)
101 icase = nint(uparam(9))
102 DO j=1,nload
103 rate(j) = uparam(9 + 2*j - 1 )
104 yfac(j) = uparam(9 + 2*j )
105 ENDDO
106 yfac_unl = uparam(9 + 2*nload + 2 )
107C initialization
108 IF(time == zero)THEN
109 uvar(1:nel,1:nuvar) = zero
110 uvar(1:nel,1) = one
111 uvar(1:nel,18)= one
112 ENDIF
113 g_ini = three_half*rbulk*(one - two*nu)/(one + nu)
114C
115C principal stretch (def gradient eigenvalues)
116 DO i=1,nel
117 trav(i) = epsxx(i)+epsyy(i)
118 rootv(i) = sqrt((epsxx(i)-epsyy(i))*(epsxx(i)-epsyy(i))
119 . + epsxy(i)*epsxy(i))
120 evv(i,1) = half*(trav(i)+rootv(i))
121 evv(i,2) = half*(trav(i)-rootv(i))
122 evv(i,3) = zero
123 ENDDO
124C rot matrix (eigenvectors)
125 DO i=1,nel
126 IF(abs(evv(i,2)-evv(i,1))<em10) THEN
127 eigv(i,1,1) = one
128 eigv(i,2,1) = one
129 eigv(i,3,1) = zero
130 eigv(i,1,2) = zero
131 eigv(i,2,2) = zero
132 eigv(i,3,2) = zero
133 ELSE
134 eigv(i,1,1) = (epsxx(i)-evv(i,2)) /rootv(i)
135 eigv(i,2,1) = (epsyy(i)-evv(i,2)) /rootv(i)
136 eigv(i,1,2) = (evv(i,1)-epsxx(i)) /rootv(i)
137 eigv(i,2,2) = (evv(i,1)-epsyy(i)) /rootv(i)
138 eigv(i,3,1) = (half*epsxy(i)) /rootv(i)
139 eigv(i,3,2) =-(half*epsxy(i)) /rootv(i)
140 ENDIF
141 ENDDO
142C Strain definition
143 IF (ismstr == 1 .OR. ismstr == 3 .OR. ismstr == 11) THEN ! engineering strain
144 DO i=1,nel
145 ev(i,1)=evv(i,1)+ one
146 ev(i,2)=evv(i,2)+ one
147 ev(i,3)=uvar(i,1)
148 ENDDO
149 ELSEIF(ismstr == 10) THEN
150 DO i=1,nel
151 ev(i,1)=sqrt(evv(i,1)+ one)
152 ev(i,2)=sqrt(evv(i,2)+ one)
153 ev(i,3)=one/ev(i,1)/ev(i,2)
154 ENDDO
155 ELSE ! true strain
156 DO i=1,nel
157 ev(i,1)=exp(evv(i,1))
158 ev(i,2)=exp(evv(i,2))
159 ev(i,3)=uvar(i,1)
160 ENDDO
161 ENDIF
162 ee(1:nel, 1) = ev(1:nel, 1) - one
163 ee(1:nel, 2) = ev(1:nel, 2) - one
164 ee(1:nel, 3) = ev(1:nel, 3) - one
165 eps_eq(1:nel) = sqrt(ee(1:nel,1)**2 + ee(1:nel,2)**2 + ee(1:nel,3)**2 )
166 deps_eq(1:nel)= eps_eq(1:nel) - uvar(1:nel,2)
167 uvar(1:nel,2)= eps_eq(1:nel)
168 gmax(1:nel) = g_ini
169 DO i=1,nel
170 f1(i) = ev(i,1)*yfac(1)*finter(ifunc(1),ee(i,1),npf,tf,df1)
171 f2(i) = ev(i,2)*yfac(1)*finter(ifunc(1),ee(i,2),npf,tf,df2)
172 !!
173 gmax(i) = max(gmax(i), yfac(1)*df1, yfac(1)*df2 )
174 fac = -half
175 dd1 = ev(i,1)**fac - one
176 nn = 1
177 DO WHILE (abs(dd1) >= em03 .AND. nn <= 20)
178 f1(i) = f1(i) +
179 . (dd1 + one)*yfac(1)*finter(ifunc(1),dd1,npf,tf,df)
180 fac = fac*(-half)
181 dd1 = ev(i,1)**fac - one
182 nn = nn + 1
183 ENDDO
184 fac = -half
185 dd2 = ev(i,2)**fac - one
186 nn = 1
187 DO WHILE (abs(dd2) >= em03 .AND. nn <= 20 )
188 f2(i) = f2(i) +
189 . (dd2 + one)*yfac(1)*finter(ifunc(1),dd2,npf,tf,df)
190 fac = fac*(-half)
191 dd2 = ev(i,2)**fac - one
192 nn = nn + 1
193 ENDDO
194 ENDDO
195!
196C--------------------------------------
197C Newton method => Find EV(3) : T3(EV(3)) = 0
198C--------------------------------------
199 DO iter = 1,5
200 DO i=1,nel
201 rv(i) = ev(i,1)*ev(i,2)*ev(i,3)
202 ee(i,3) = ev(i,3) - one
203 y3 = yfac(1)*finter(ifunc(1),ee(i,3),npf,tf,df3)
204 f3(i) = ev(i,3)*y3
205 gmax(i) = max(gmax(i), yfac(1)*df3 )
206 dflam3 = zero
207 fac = -half
208 scale = ev(i,3)**fac
209 dd3 = scale - one
210 nn = 1
211 DO WHILE (abs(dd3) >= em03 .AND. nn <= 20 )
212 y33 = yfac(1)*finter(ifunc(1),dd3,npf,tf,df)
213 f3(i) = f3(i) + scale*y33
214 dflam3 = dflam3 + y33*fac * scale/ev(i,3) + yfac(1)*fac*scale**2*df/ev(i,3)
215 fac = fac*(-half)
216 scale = ev(i,3)**fac
217 dd3 = scale - one
218 nn = nn + 1
219 END DO
220
221 invr = one/max(em20,rv(i))
222 p = invr*rbulk*(rv(i) - one)
223 t1(i) = invr*(two_third*f1(i) - third*(f2(i) + f3(i))) + p
224 t2(i) = invr*(two_third*f2(i) - third*(f1(i) + f3(i))) + p
225 t3(i) = invr*(two_third*f3(i) - third*(f1(i) + f2(i))) + p
226C
227 kt = y3 + ev(i,3)*yfac(1)*df3 + dflam3
228 kt = two_third*kt*invr + rbulk*invr/ev(i,3) -
229 . invr*(two_third*f3(i) - third*(f1(i)+f2(i)))/ev(i,3)
230 IF(kt /= zero) ev(i,3) = ev(i,3) - t3(i)/kt
231 ENDDO ! NEL
232 ENDDO ! ITER
233 ne_l = 0
234 ne_unl = 0
235 DO i=1,nel
236 ecurent(i)= half*(ee(i,1)*t1(i) + ee(i,2)*t2(i))
237 uvar(i,17)= max(uvar(i,17), ecurent(i))
238 deint(i) = ecurent(i) - uvar(i,16)
239 uvar(i,16) = ecurent(i)
240 IF(deint(i) >= zero .OR. deps_eq(i) >= zero) THEN
241 ne_l = ne_l + 1
242 indx_l(ne_l) = i
243 uvar(i,19) = 1
244 IF(uvar(i,18) == one )uvar(i,17) = uvar(i,16)
245 ELSE
246 ne_unl = ne_unl + 1
247 indx_unl(ne_unl) = i
248 ENDIF
249 ENDDO
250 DO i=1,nel
251 IF(uvar(i,17) > zero .AND. ecurent(i) <= uvar(i,17)) THEN
252 dam = one - (ecurent(i)/uvar(i,17))**shape
253 dam = one - (one - hys)*dam
254C global
255 t1(i) = dam*t1(i)
256 t2(i) = dam*t2(i)
257 uvar(i,19) = -1
258 uvar(i,18) = dam
259 ENDIF
260 ENDDO
261C-------------------------------------------------------------
262C-------------------------------------------------------------
263 IF (ismstr == 1 .OR. ismstr == 3 .OR. ismstr == 11) THEN ! engineering strain
264 DO i=1,nel
265 epszz(i) =ev(i,3) - one
266 uvar(i,1) = ev(i,3)
267 ENDDO
268 ELSEIF (ismstr == 10) THEN ! left gauchy-green strain
269 DO i=1,nel
270 epszz(i) =ev(i,3) - one
271 uvar(i,1) = ev(i,3)
272 ENDDO
273 ELSE ! true strain
274 DO i=1,nel
275 epszz(i) =log(ev(i,3))
276 uvar(i,1) = ev(i,3)
277 ENDDO
278 ENDIF
279 DO i=1,nel
280 rv(i) = ev(i,1)*ev(i,2)*ev(i,3)
281 dezz(i) =-nu/(one-nu)*(depsxx(i)+depsyy(i))
282 signxx(i) =eigv(i,1,1)*t1(i)+eigv(i,1,2)*t2(i)
283 signyy(i) =eigv(i,2,1)*t1(i)+eigv(i,2,2)*t2(i)
284 signxy(i) =eigv(i,3,1)*t1(i)+eigv(i,3,2)*t2(i)
285C
286 signyz(i) = sigoyz(i)+gs(i)*depsyz(i)
287 signzx(i) = sigozx(i)+gs(i)*depszx(i)
288 thkn(i) = thkn(i) + dezz(i)*thklyl(i)*off(i)
289 viscmax(i)= zero
290 emax = two*g_ini*(one + nu)
291 a11 = emax/(one - nu**2)
292 soundsp(i)= sqrt(a11/rho0(i))
293 ENDDO
294C-----------
295 RETURN
296 END
max
#define max(a, b)
Definition macros.h:21
sigeps88c
subroutine sigeps88c(nel, nuparam, nuvar, nfunc, ifunc, npf, npt0, ilayer, ngl, off, ismstr, gs, tf, time, timestep, uparam, rho0, area, eint, thklyl, depsxx, depsyy, depsxy, depsyz, depszx, epsxx, epsyy, epsxy, epsyz, epszx, sigoxx, sigoyy, sigoxy, sigoyz, sigozx, signxx, signyy, signxy, signyz, signzx, soundsp, viscmax, thkn, uvar)
Definition sigeps88c.F:40