OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps111.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!|| sigeps111 ../engine/source/materials/mat/mat111/sigeps111.F
25!||--- called by ------------------------------------------------------
26!|| mulaw ../engine/source/materials/mat_share/mulaw.F90
27!||--- calls -----------------------------------------------------
28!|| cardan_method ../engine/source/materials/mat/mat111/sigeps111.F
29!|| finter ../engine/source/tools/curve/finter.F
30!|| prodaat ../engine/source/materials/tools/prodAAT.F
31!||====================================================================
32 SUBROUTINE sigeps111(
33 1 NEL , NUPARAM, NUVAR , NFUNC , IFUNC , NPF ,
34 2 TF , TIME , TIMESTEP, UPARAM, RHO0 , RHO ,
35 3 VOLUME , EINT , NGL ,
36 4 EPSPXX , EPSPYY , EPSPZZ , EPSPXY, EPSPYZ, EPSPZX,
37 5 DEPSXX , DEPSYY , DEPSZZ , DEPSXY, DEPSYZ, DEPSZX,
38 6 EPSXX , EPSYY , EPSZZ , EPSXY , EPSYZ , EPSZX ,
39 7 SIGOXX , SIGOYY , SIGOZZ , SIGOXY, SIGOYZ, SIGOZX,
40 8 SIGNXX , SIGNYY , SIGNZZ , SIGNXY, SIGNYZ, SIGNZX,
41 9 SIGVXX , SIGVYY , SIGVZZ , SIGVXY, SIGVYZ, SIGVZX,
42 A SOUNDSP, VISCMAX, UVAR , OFF , ISMSTR, ET ,
43 B IHET , OFFG , EPSTH3 , IEXPAN ,MFXX ,MFXY ,
44 C MFXZ , MFYX , MFYY , MFYZ ,MFZX ,MFZY ,
45 D MFZZ )
46C-----------------------------------------------
47C I M P L I C I T T Y P E S
48C-----------------------------------------------
49#include "implicit_f.inc"
50C----------------------------------------------------------------
51C I N P U T A R G U M E N T S
52C----------------------------------------------------------------
53 INTEGER :: NEL, NUPARAM, NUVAR,ISMSTR,NGL(*),IHET,IEXPAN
54 my_real :: TIME , TIMESTEP , UPARAM(NUPARAM)
55 my_real , DIMENSION(NEL) :: RHO , RHO0 , VOLUME, EINT,
56 . EPSPXX, EPSPYY, EPSPZZ, EPSPXY,
57 . EPSPYZ, EPSPZX, DEPSXX, DEPSYY,
58 . DEPSZZ, DEPSXY, DEPSYZ, DEPSZX,
59 . EPSXX , EPSYY , EPSZZ , EPSXY ,
60 . EPSYZ , EPSZX , SIGOXX, SIGOYY,
61 . SIGOZZ, SIGOXY, SIGOYZ, SIGOZX,
62 . OFFG , EPSTH3,
63 . MFXX , MFXY , MFXZ , MFYX ,
64 . MFYY , MFYZ , MFZX , MFZY ,
65 . MFZZ
66C----------------------------------------------------------------
67C O U T P U T A R G U M E N T S
68C----------------------------------------------------------------
69 my_real , DIMENSION(NEL) :: signxx , signyy , signzz,
70 . signxy , signyz , signzx,
71 . sigvxx , sigvyy , sigvzz,
72 . sigvxy , sigvyz , sigvzx,
73 . soundsp, viscmax, et
74C----------------------------------------------------------------
75C I N P U T O U T P U T A R G U M E N T S
76C----------------------------------------------------------------
77 my_real , DIMENSION(NEL) :: off
78 my_real , DIMENSION(NEL,NUVAR) :: uvar
79C----------------------------------------------------------------
80C VARIABLES FOR FUNCTION INTERPOLATION
81C----------------------------------------------------------------
82 INTEGER :: NPF(*), NFUNC, IFUNC(NFUNC)
83 my_real :: FINTER,FINTTE,TF(*),FINT2V
84 EXTERNAL FINTER,FINTTE
85C----------------------------------------------------------------
86C L O C A L V A R I B L E S
87C----------------------------------------------------------------
88 INTEGER :: I,J,ITEST,ONLY_TENSION_DATA
89 INTEGER , DIMENSION(NEL) :: ICOMP
90 my_real :: dw_di,l2,l3,l5,delta,
91 . g,rbulk,scale,t,epst,df,du,
92 . trace,p,fac,invjdetf,nu,gini,rbulkini,diff
93 my_real :: b(nel,6),lam2(nel,3),
94 . lambda(nel,3),evd(3),bd(nel,6)
95 my_real ,DIMENSION(NEL) :: jdetf,i1,i1b,lamt,j2third,
96 . t1,t2,t3,eti,gs,k
97 my_real, DIMENSION(NEL,3,3) :: f,fft
98c----------------------------------------------------------------
99c UVAR(1) = engineering strain
100c UVAR(2) = engineering stress ! for unixial test only for now
101c----------------------------------------------------------------
102c material parameters
103 nu = uparam(1)
104 itest = nint(uparam(2))
105 scale = uparam(3)
106 g = uparam(4)
107 rbulk = uparam(5)
108 gini = uparam(6)
109 only_tension_data = nint(uparam(7)) ! == 1 existing only tension data
110C
111 DO i=1,nel
112 f(i,1,1) = one + mfxx(i)
113 f(i,2,2) = one + mfyy(i)
114 f(i,3,3) = one + mfzz(i)
115 f(i,1,2) = mfxy(i)
116 f(i,2,3) = mfyz(i)
117 f(i,3,1) = mfzx(i)
118 f(i,2,1) = mfyx(i)
119 f(i,3,2) = mfzy(i)
120 f(i,1,3) = mfxz(i)
121 ENDDO
122C
123 CALL prodaat(f , fft, nel) ! B = F * FT
124
125 DO i=1,nel
126 b(i,1) = fft(i,1,1)
127 b(i,2) = fft(i,2,2)
128 b(i,3) = fft(i,3,3)
129 b(i,4) = fft(i,1,2)
130 b(i,5) = fft(i,2,3)
131 b(i,6) = fft(i,3,1)
132 ENDDO
133 icomp(1:nel) = 0
134 DO i = 1,nel
135 ! JDETF = RHO0/RHO = RELATIVE VOLUME = DET F (F = GRADIENT OF DEFORMATION)
136 jdetf(i) = f(i,1,1)*f(i,2,2)*f(i,3,3) - f(i,1,1)*f(i,2,3)*f(i,3,2)
137 . - f(i,3,3)*f(i,1,2)*f(i,2,1) + f(i,1,2)*f(i,2,3)*f(i,3,1)
138 . + f(i,2,1)*f(i,3,2)*f(i,1,3) - f(i,2,2)*f(i,3,1)*f(i,1,3)
139 !FIRST INVARIANT of B=FFT
140 i1b(i) = fft(i,1,1) + fft(i,2,2) + fft(i,3,3)
141C
142 IF (jdetf(i) > zero) THEN
143 j2third(i) = exp((-two_third )*log(jdetf(i)))
144 ELSE
145 j2third(i) = zero
146 ENDIF
147C
148 trace = (b(i,1)+b(i,2)+ b(i,3))*third
149 bd(i,1) = j2third(i)*(b(i,1) - trace)
150 bd(i,2) = j2third(i)*(b(i,2) - trace)
151 bd(i,3) = j2third(i)*(b(i,3) - trace)
152 bd(i,4:6) = j2third(i)*b(i,4:6)
153 i1(i) = j2third(i)*i1b(i)
154 ENDDO
155
156 IF(only_tension_data == 0) THEN
157 DO i = 1,nel
158 diff = jdetf(i) - one
159 IF(diff < zero) icomp(i) = 1
160 ENDDO
161 ENDIF
162C
163 SELECT CASE (itest)
164
165 CASE (1) ! Uniaxial traction test (tension)
166 !
167 ! solving lt**3 - I*lt +2 = 0
168 !
169 !!CALL CARDAN_METHOD(NEL,I1B,LAMT)
170 CALL cardan_method(itest,nel,i1,icomp,lamt)
171
172 DO i=1,nel
173 invjdetf = one/max(em20, jdetf(i))
174 epst = lamt(i) - one
175 t = scale*finter(ifunc(1),epst,npf,tf,df)
176 l2 = lamt(i)**2
177 dw_di = zero
178 IF (i1b(i) == three .OR. lamt(i) == one) THEN
179 dw_di = zero
180 ELSEIF((lamt(i)*l2 - one) /= zero)THEN
181 dw_di = t*l2 / (lamt(i)*l2 - one)
182 ENDIF
183 gs(i)= max(gini,df*scale)
184 k(i) = two*gs(i)*(one+nu) / max(em20,three*(one-two*nu))
185 p = rbulk*(jdetf(i) - one)
186 signxx(i) = invjdetf*dw_di * bd(i,1) + p
187 signyy(i) = invjdetf*dw_di * bd(i,2) + p
188 signzz(i) = invjdetf*dw_di * bd(i,3) + p
189 signxy(i) = invjdetf*dw_di * bd(i,4)
190 signyz(i) = invjdetf*dw_di * bd(i,5)
191 signzx(i) = invjdetf*dw_di * bd(i,6)
192 uvar(i,1) = lamt(i) - one
193 uvar(i,2) = jdetf(i)*signzz(i) / lamt(i) ! only for testing
194 ENDDO
195 CASE(2) ! equibiaxail traction test
196 !
197 ! solving 2*lt**6 - I*lt**4 + 1 = 0
198 ! By considering lt**2 = x --> 2*x**3 - I*x**2 +1 = 0
199 ! by considering z = 1/x we get z**3 - I*Z + 2 = 0
200 ! Z ---> lt = sqrt(1/z)
201!!
202 CALL cardan_method(itest,nel,i1,icomp,lamt)
203!!
204 DO i=1,nel
205 invjdetf = one/max(em20, jdetf(i))
206 lamt(i) = one/max(em20,lamt(i))
207 lamt(i) = sqrt(lamt(i))
208 epst = lamt(i) - one
209 t = scale*finter(ifunc(1),epst,npf,tf,df)
210 l5 = lamt(i)**(-5)
211 dw_di = zero
212 IF (i1b(i) == three .OR. lamt(i) == one) THEN
213 dw_di = zero
214 ELSEIF((lamt(i) - l5) /= zero) THEN
215 dw_di = t / (lamt(i) - l5)
216 ENDIF
217 gs(i)= max(gini,df*scale)
218 k(i) = two*gs(i)*(one+nu) / max(em20,three*(one-two*nu))
219 p = rbulk*(jdetf(i) - one)
220 signxx(i) = invjdetf*dw_di * bd(i,1) + p
221 signyy(i) = invjdetf*dw_di * bd(i,2) + p
222 signzz(i) = invjdetf*dw_di * bd(i,3) + p
223 signxy(i) = invjdetf*dw_di * bd(i,4)
224 signyz(i) = invjdetf*dw_di * bd(i,5)
225 signzx(i) = invjdetf*dw_di * bd(i,6)
226 uvar(i,1) = lamt(i) - one
227 uvar(i,2) = jdetf(i)*signzz(i) / lamt(i)
228 ENDDO
229 CASE(3) ! planar traction test
230
231 !
232 ! solving lt**4 +(1- I)*lt**2 + 1 = 0
233 ! By considering lt**2 = x --> x**2 +(1-I)*x +1 = 0
234 ! lt = sqrt(X)
235!!
236 DO i=1,nel
237 delta = (one - i1(i)) ** 2 - four
238 lamt(i) = one
239 IF(delta > 0) THEN
240 lamt(i) = half*max(i1(i)+ sqrt(delta) - one ,i1(i) - sqrt(delta) - one)
241 ELSEIF(delta == zero) THEN
242 lamt(i) = half*(i1(i) - one)
243 ENDIF
244 ENDDO
245!!
246 DO i=1,nel
247 invjdetf= one/max(em20, jdetf(i))
248 lamt(i) = sqrt(max(zero,lamt(i)))
249 epst = lamt(i) - one
250 t = scale*finter(ifunc(1),epst,npf,tf,df)
251 l3 = lamt(i)**(-3)
252 dw_di = zero
253 IF (i1b(i) == three .OR. lamt(i) == one) THEN
254 dw_di = zero
255 ELSEIF((lamt(i) - l3) /= zero) THEN
256 dw_di = t / (lamt(i) - l3)
257 ENDIF
258 gs(i)= max(gini,df*scale)
259 k(i) = two*gs(i)*(one+nu) / max(em20,three*(one-two*nu))
260 p = rbulk*(jdetf(i) - one)
261 signxx(i) = invjdetf*dw_di * bd(i,1) + p
262 signyy(i) = invjdetf*dw_di * bd(i,2) + p
263 signzz(i) = invjdetf*dw_di * bd(i,3) + p
264 signxy(i) = invjdetf*dw_di * bd(i,4)
265 signyz(i) = invjdetf*dw_di * bd(i,5)
266 signzx(i) = invjdetf*dw_di * bd(i,6)
267 uvar(i,1) = lamt(i) - one
268 uvar(i,2) = -jdetf(i)*signzx(i) / lamt(i) ! only for testing
269 ENDDO
270 END SELECT
271C
272 IF (ihet > 1) THEN
273 DO i=1,nel
274 et(i) = max(one, gs(i)/gini)
275 ENDDO
276 ENDIF
277C----------------
278 DO i=1,nel
279 soundsp(i) = sqrt((four_over_3*g + rbulk)/rho(i))
280 viscmax(i) = zero
281 ENDDO
282c-----------
283 RETURN
284 END
285C
286!||====================================================================
287!|| cardan_method ../engine/source/materials/mat/mat111/sigeps111.F
288!||--- called by ------------------------------------------------------
289!|| sigeps111 ../engine/source/materials/mat/mat111/sigeps111.F
290!||====================================================================
291 SUBROUTINE cardan_method(ITEST,NEL,I1,ICOMP,LAMT)
292C solve X**3 - I*X + 2 = 0
293C-----------------------------------------------
294C I m p l i c i t T y p e s
295C-----------------------------------------------
296#include "implicit_f.inc"
297C-----------------------------------------------
298C G l o b a l V a r i a b l e s
299C-----------------------------------------------
300 INTEGER :: ITEST,NEL
301 INTEGER :: ICOMP(NEL)
302 my_real , DIMENSION(NEL):: I1(NEL),LAMT(NEL)
303C-----------------------------------------------
304C L o c a l V a r i a b l e s
305C-----------------------------------------------
306 INTEGER :: I
307 my_real
308 . X0,X1,X2,DELTA,AA,BB,HUNDRED8,HUIT1,ROOT,FACT
309C
310 HUNDRED8 = hundred + eight
311 huit1 = 81
312 IF(itest == 1) THEN
313 DO i=1,nel
314 delta = four*i1(i)**3 - hundred8
315 IF(delta > zero) THEN
316 aa = three/i1(i)
317 bb = third*i1(i)
318 aa = -aa*sqrt(aa)
319 x0 = two*sqrt(bb)*cos(third*acos(aa))
320 x1 = two*sqrt(bb)*cos(third*acos(aa) + two_third*pi)
321 x2 = two*sqrt(bb)*cos(third*acos(aa) + four_over_3*pi)
322 lamt(i) = max(one,x0,x1,x2)! tension
323 IF(icomp(i) == 0) cycle
324 ! in compression
325 IF(x0*x1 < zero) THEN
326 IF(x0 > zero) lamt(i) = x0
327 IF(x1 > zero) lamt(i) = x1
328 IF(lamt(i) > x2 .AND. x2 > zero ) lamt(i) = x2
329 ELSE
330 IF(x0 < zero) THEN
331 IF(x2 > zero) lamt(i) = x2
332 ELSE
333 lamt(i) = min(x0, x1)
334 IF(lamt(i) > x2 .AND. x2 > zero ) lamt(i) = x2
335 ENDIF
336 ENDIF
337 ELSE
338 lamt(i) = max(one,three/i1(i),-six/i1(i))
339 ENDIF
340 ENDDO
341 ELSEIF(itest == 2) THEN
342 DO i=1,nel
343 delta = four*i1(i)**3 - hundred8
344 IF(delta > zero) THEN
345 aa = three/i1(i)
346 bb = third*i1(i)
347 aa = -aa*sqrt(aa)
348 x0 = two*sqrt(bb)*cos(third*acos(aa))
349 x1 = two*sqrt(bb)*cos(third*acos(aa) + two_third*pi)
350 x2 = two*sqrt(bb)*cos(third*acos(aa) + four_over_3*pi)
351 lamt(i) = max(x0,x1,x2) ! COMPRESSION
352 IF(icomp(i) == 1) cycle
353 ! In tension
354 IF(x0*x1 < zero) THEN
355 IF(x0 > zero) lamt(i) = x0
356 IF(x1 > zero) lamt(i) = x1
357 IF(lamt(i) > x2 .AND. x2 > zero ) lamt(i) = x2
358 ELSE
359 IF(x0 < zero) THEN
360 IF(x2 > zero) lamt(i) = x2
361 ELSE
362 lamt(i) = min(x0, x1)
363 IF(lamt(i) > x2 .AND. x2 > zero ) lamt(i) = x2
364 ENDIF
365 ENDIF
366 ELSE
367 lamt(i) = max(three/i1(i),-six/i1(i))
368 ENDIF
369c
370 ENDDO
371
372
373
374 ENDIF
375C-----
376 RETURN
377 END SUBROUTINE cardan_method
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
prodaat
subroutine prodaat(a, c, nel)
Definition prodAAT.F:34
sigeps111
subroutine sigeps111(nel, nuparam, nuvar, nfunc, ifunc, npf, tf, time, timestep, uparam, rho0, rho, volume, eint, ngl, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, depsxx, depsyy, depszz, depsxy, depsyz, depszx, epsxx, epsyy, epszz, epsxy, epsyz, epszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, sigvxx, sigvyy, sigvzz, sigvxy, sigvyz, sigvzx, soundsp, viscmax, uvar, off, ismstr, et, ihet, offg, epsth3, iexpan, mfxx, mfxy, mfxz, mfyx, mfyy, mfyz, mfzx, mfzy, mfzz)
Definition sigeps111.F:46
cardan_method
subroutine cardan_method(itest, nel, i1, icomp, lamt)
Definition sigeps111.F:292