OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
law119_membrane.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!|| law119_membrane ../engine/source/materials/mat/mat119/law119_membrane.F
25!||--- called by ------------------------------------------------------
26!|| sigeps119c ../engine/source/materials/mat/mat119/sigeps119c.F
27!||--- calls -----------------------------------------------------
28!|| table_vinterp ../engine/source/tools/curve/table_tools.F
29!||--- uses -----------------------------------------------------
30!|| interface_table_mod ../engine/share/modules/table_mod.F
31!|| message_mod ../engine/share/message_module/message_mod.F
32!|| table_mod ../engine/share/modules/table_mod.F
33!||====================================================================
34 SUBROUTINE law119_membrane(
35 . NEL ,NUPARAM,NUVAR ,UPARAM ,UVAR ,
36 . GS ,ET ,DEPSXX ,DEPSYY ,DEPSXY ,
37 . EPSXX ,EPSYY ,EPSXY ,EPSYZ ,EPSZX ,
38 . SIGOXX ,SIGOYY ,SIGOXY ,SIGOYZ ,SIGOZX ,
39 . SIGNXX ,SIGNYY ,SIGNXY ,SIGNYZ ,SIGNZX ,
40 . NUMTABL,ITABLE ,TABLE )
41C-----------------------------------------------
42C M o d u l e s
43C-----------------------------------------------
44 USE table_mod
45 USE message_mod
46 USE interface_table_mod
47C-----------------------------------------------
48C I m p l i c i t T y p e s
49C-----------------------------------------------
50#include "implicit_f.inc"
51C-----------------------------------------------
52C C o m m o n B l o c k s
53C-----------------------------------------------
54#include "com04_c.inc"
55C-----------------------------------------------
56C D u m m y A r g u m e n t s
57C-----------------------------------------------
58 INTEGER ,INTENT(IN) :: NEL,NUPARAM,NUVAR,NUMTABL
59 INTEGER ,DIMENSION(NUMTABL),INTENT(IN) :: ITABLE
60 my_real ,INTENT(IN) :: GS
61 my_real ,DIMENSION(NUPARAM) ,INTENT(IN) :: UPARAM
62 my_real ,DIMENSION(NEL) ,INTENT(IN) :: DEPSXX,DEPSYY,DEPSXY,
63 . EPSXX,EPSYY,EPSXY,EPSYZ,EPSZX,SIGOXX,SIGOYY,SIGOXY,SIGOYZ,SIGOZX
64 my_real ,DIMENSION(NEL) ,INTENT(OUT) ::
65 . signxx,signyy,signxy,signyz,signzx,et
66 my_real ,DIMENSION(NEL,NUVAR) ,INTENT(INOUT) :: uvar
67 TYPE(ttable), DIMENSION(NTABLE) :: TABLE
68C-----------------------------------------------
69C L o c a l V a r i a b l e s
70C-----------------------------------------------
71 INTEGER :: I,II,NINDXU,NINDXL,NINDXR,NINDXT,NINDXC,NREAC,FUNC,FUND,FUNCR,
72 . ndim,ireload
73 my_real :: nu12,nu21,g12,a11,a12,a22,rcomp,
74 . fscale,fscale1,fscale2,fscalet,dw,ht,hf,
75 . det,s,d,t,p1,p2,r,dezz,svm2,xint,yint
76 INTEGER ,DIMENSION(NEL) :: INDXU,INDXL,INDXR,INDXT,INDXC,INDXX
77 my_real ,DIMENSION(NEL) :: epsq,svm,etx,etl,etu,exx,eyy,dx,dy,
78 . sigf,emax,smax,eminrl,sminrl,emaxrl,smaxrl,beta
79 INTEGER ,DIMENSION(NEL,1) :: IPOS1
80 my_real ,DIMENSION(NEL,1) :: XX1
81C-----------------------------------------------
82C UVAR(1) = Emax
83C UVAR(2) = Smax
84C UVAR(3) = EminRL
85C UVAR(4) = EmaxRL
86C UVAR(5) = SminRL
87C UVAR(6) = SmaxRL
88C UVAR(7) = reactivation flag
89C UVAR(8) = EPSQ - equivalent strain
90C UVAR(9) =
91C UVAR(10)= compression flag
92C=======================================================================
93 nu12 = uparam(3)
94 nu21 = uparam(4)
95 g12 = uparam(6)
96 a11 = uparam(7)
97 a22 = uparam(8)
98 a12 = uparam(9)
99 rcomp = uparam(10)
100 fscale1 = uparam(11)
101 fscale2 = uparam(12)
102 fscalet = uparam(13)
103 xint = uparam(18)
104 yint = uparam(19)
105 ireload = nint(uparam(21))
106 det = one / (one - nu12*nu21)
107c
108 func = itable(1)
109 fund = itable(2)
110 nindxt = 0
111 nindxc = 0
112 nindxu = 0
113 nindxl = 0
114 nindxr = 0
115 nreac = 0
116c-----------------------------------------
117 emax(1:nel) = uvar(1:nel,1)
118 smax(1:nel) = uvar(1:nel,2)
119 eminrl(1:nel) = uvar(1:nel,3)
120 emaxrl(1:nel) = uvar(1:nel,4)
121 sminrl(1:nel) = uvar(1:nel,5)
122 smaxrl(1:nel) = uvar(1:nel,6)
123c-----------------------------------------
124c test principal strain direction : tag elements in tension / compression
125c-----------------------------------------
126 DO i=1,nel
127 s = half*(epsxx(i) + epsyy(i))
128 d = half*(epsxx(i) - epsyy(i))
129 r = sqrt(epsxy(i)**2 + d*d)
130 p1 = s + r
131 p2 = s - r
132 IF (p1 > zero .and. p1 >= -p2) THEN ! max principal deformation in tension
133 IF (nint(uvar(i,7)) == 1) THEN ! elements are reactivated after passing through slipring
134 nreac = nreac + 1
135 indxx(nreac) = i
136 ELSE
137 nindxt = nindxt + 1
138 indxt(nindxt) = i
139 END IF
140 beta(i) = one
141 ELSE
142 nindxc = nindxc + 1
143 indxc(nindxc) = i
144 beta(i) = rcomp
145 END IF
146 END DO
147c-----------------------------------------
148c Shear and compression - always linear using total strain
149c-----------------------------------------
150 DO i=1,nel
151 signxy(i) = g12*epsxy(i) * beta(i)
152 signyz(i) = gs *epsyz(i) * beta(i)
153 signzx(i) = gs *epszx(i) * beta(i)
154 et(i) = beta(i)
155 END DO
156c-----------------------------------------
157c In plane stress : linear using total strain when func = 0 or compression
158c-----------------------------------------
159 IF (func == 0) THEN
160 DO i=1,nel
161 signxx(i) = (a11*epsxx(i) + a12*epsyy(i))*beta(i)
162 signyy(i) = (a12*epsxx(i) + a22*epsyy(i))*beta(i)
163 END DO
164c
165 ELSE ! FUNC > 0
166c
167 IF (nindxt > 0 .and. fund == 0) THEN ! nonlinear tension
168 ndim = table(func)%NDIM
169 DO i=1,nel
170 epsq(i) = sqrt((epsxx(i)**2 + epsyy(i)**2) / (one + nu21**2))
171 ENDDO
172 xx1(1:nel,1) = epsq(1:nel)
173 ipos1(1:nel,1) = 1
174c
175 CALL table_vinterp(table(func),nel,nel,ipos1,xx1,sigf,etl)
176c
177 DO ii=1,nindxt
178 i = indxt(ii)
179 a11 = etl(i) * det * fscale1
180 a22 = a11 * fscalet
181 a12 = a11 * nu21
182 signxx(i) = sigoxx(i) + a11*depsxx(i) + a12*depsyy(i)
183 signyy(i) = sigoyy(i) + a12*depsxx(i) + a22*depsyy(i)
184 ENDDO
185 END IF
186c
187 IF (nindxc > 0) THEN ! compression
188 DO ii=1,nindxc
189 i = indxc(ii)
190 signxx(i) = (a11*epsxx(i) + a12*epsyy(i))*rcomp
191 signyy(i) = (a12*epsxx(i) + a22*epsyy(i))*rcomp
192 END DO
193 IF (fund > 0) THEN
194 DO ii=1,nindxc
195 i = indxc(ii)
196 emax(i) = em20
197 smax(i) = zero
198 eminrl(i) = zero
199 emaxrl(i) = em20
200 sminrl(i) = zero
201 smaxrl(i) = zero
202 END DO
203 END IF
204 END IF
205c
206 IF (nreac > 0) THEN ! reactivated elements after slipring
207 DO i=1,nel
208 epsq(i) = sqrt((epsxx(i)**2 + epsyy(i)**2) / (one + nu21**2))
209 ENDDO
210 xx1(1:nel,1) = epsq(1:nel)
211 ipos1(1:nel,1) = 1
212c
213 CALL table_vinterp(table(func),nel,nel,ipos1,xx1,sigf,etl)
214c
215 DO ii=1,nreac
216 i = indxx(ii)
217 signxx(i) = (a11*epsxx(i) + a12*epsyy(i))*beta(i)
218 signyy(i) = (a12*epsxx(i) + a22*epsyy(i))*beta(i)
219 svm2 = signxx(i)**2 + signyy(i)**2 - signxx(i)*signyy(i)
220 svm(i) = sqrt(svm2)
221 emax(i) = epsq(i)
222 smax(i) = svm(i)
223 eminrl(i) = zero
224 emaxrl(i) = zero
225 sminrl(i) = zero
226 smaxrl(i) = zero
227 END DO
228 END IF
229 END IF
230c-----------------------------------------
231c
232 IF (func > 0) THEN
233 ndim = table(func)%NDIM
234 DO i=1,nel
235 epsq(i) = sqrt((epsxx(i)**2 + epsyy(i)**2) / (one + nu21**2))
236 ENDDO
237c
238 xx1(1:nel,1) = epsq(1:nel)
239 ipos1(1:nel,1) = 1
240c
241 CALL table_vinterp(table(func),nel,nel,ipos1,xx1,sigf,etl)
242c
243 DO ii=1,nindxt
244 i = indxt(ii)
245 a11 = etl(i) * det * fscale1
246 a22 = a11 * fscalet
247 a12 = a11 * nu21
248 signxx(i) = sigoxx(i) + a11*depsxx(i) + a12*depsyy(i)
249 signyy(i) = sigoyy(i) + a12*depsxx(i) + a22*depsyy(i)
250 ENDDO
251 END IF
252c-----------------------------------------
253c
254 IF (fund > 0) THEN ! nonlinear loading and unloading with hysteresis
255c
256 DO i=1,nel
257 epsq(i) = sqrt((epsxx(i)**2 + epsyy(i)**2) / (one + nu21**2))
258 svm(i) = sqrt(signxx(i)**2 + signyy(i)**2 - signxx(i)*signyy(i))
259 ENDDO
260c
261 DO ii=1,nindxt
262 i = indxt(ii)
263 dw = epsq(i) - uvar(i,8)
264 IF (dw < zero .and. uvar(i,10) >= zero) THEN
265 nindxu = nindxu + 1
266 indxu(nindxu) = i
267 ELSE IF (svm(i) >= smax(i) .or. uvar(i,10) == -one) THEN ! loading
268 nindxl = nindxl + 1
269 indxl(nindxl) = i
270 ELSE ! reloading
271 nindxr = nindxr + 1
272 indxr(nindxr) = i
273 END IF
274 END DO
275c loading
276 IF (nindxl > 0) THEN
277 DO ii=1,nindxl
278 i = indxl(ii)
279 a11 = etl(i) * det * fscale1
280 a22 = a11 * fscalet
281 a12 = a11 * nu21
282 signxx(i) = sigoxx(i) + a11*depsxx(i) + a12*depsyy(i)
283 signyy(i) = sigoyy(i) + a12*depsxx(i) + a22*depsyy(i)
284 svm2 = signxx(i)**2 + signyy(i)**2 - signxx(i)*signyy(i)
285 svm(i) = sqrt(svm2)
286 emax(i) = max(em20, epsq(i))
287 eminrl(i) = epsq(i)
288 emaxrl(i) = max(em20, epsq(i))
289 smax(i) = svm(i)
290 sminrl(i) = svm(i)
291 smaxrl(i) = svm(i)
292 ENDDO
293 END IF
294c reloading
295 IF (nindxr > 0) THEN
296 IF (ireload == 0) THEN
297 funcr = fund! reloading follows unloading curve
298 fscale = fscale2
299 DO i=1,nel
300 xx1(i,1) = epsq(i) * xint / emaxrl(i)
301 END DO
302 ELSE
303 funcr = func! reloading follows loading curve
304 fscale = fscale1
305 DO i=1,nel
306 xx1(i,1) = emax(i) * (epsq(i) - eminrl(i)) / (emax(i) - eminrl(i))
307 END DO
308 END IF
309
310 ipos1(1:nel,1) = 1
311 CALL table_vinterp(table(funcr),nel,nel,ipos1,xx1,sigf,etx)
312c
313 IF (ireload == 1) then! reloading follows loading curve
314 DO ii=1,nindxr
315 i = indxr(ii)
316 ht = (smax(i)-sminrl(i)) / (emax(i)-eminrl(i))
317 hf = smax(i) / emax(i)
318 etx(i) = fscale * etx(i) * ht / hf
319c
320 a11 = etx(i) * det
321 a22 = a11 * fscalet
322 a12 = a11 * nu21
323 signxx(i) = sigoxx(i) + a11*depsxx(i) + a12*depsyy(i)
324 signyy(i) = sigoyy(i) + a12*depsxx(i) + a22*depsyy(i)
325 svm2 = signxx(i)**2 + signyy(i)**2 - signxx(i)*signyy(i)
326 svm(i) = sqrt(svm2)
327 emaxrl(i) = max(em20, epsq(i))
328 smaxrl(i) = svm(i)
329 END DO
330 ELSE ! reloading follows unloading curve
331 DO ii=1,nindxr
332 i = indxr(ii)
333 ht = smax(i) / emax(i)
334 hf = yint / xint
335 etx(i) = fscale * etx(i) * ht / hf
336c
337 a11 = etx(i) * det
338 a22 = a11 * fscalet
339 a12 = a11 * nu21
340 signxx(i) = sigoxx(i) + a11*depsxx(i) + a12*depsyy(i)
341 signyy(i) = sigoyy(i) + a12*depsxx(i) + a22*depsyy(i)
342 svm2 = signxx(i)**2 + signyy(i)**2 - signxx(i)*signyy(i)
343 svm(i) = sqrt(svm2)
344 emaxrl(i) = max(em20, emax(i))
345 smaxrl(i) = smax(i)
346 END DO
347 END IF
348 END IF ! reloading
349c unloading
350 IF (nindxu > 0) THEN
351 DO i=1,nel
352 xx1(i,1) = epsq(i) * xint / emaxrl(i)
353 END DO
354 ipos1(1:nel,1) = 1
355
356 CALL table_vinterp(table(fund),nel,nel,ipos1,xx1,sigf,etu)
357c
358 DO ii=1,nindxu
359 i = indxu(ii)
360 etx(i) = etu(i) * (smaxrl(i) / emaxrl(i)) * (xint / yint)
361 IF (epsq(i) > zero) THEN
362 etx(i) = max(etx(i), svm(i) / epsq(i))
363 END IF
364 END DO
365c
366 DO ii=1,nindxu
367 i = indxu(ii)
368 a11 = etx(i) * det * fscale2
369 a22 = a11 * fscalet
370 a12 = a11 * nu21
371 signxx(i) = sigoxx(i) + a11*depsxx(i) + a12*depsyy(i)
372 signyy(i) = sigoyy(i) + a12*depsxx(i) + a22*depsyy(i)
373 svm2 = signxx(i)**2 + signyy(i)**2 - signxx(i)*signyy(i)
374 svm(i) = sqrt(svm2)
375 eminrl(i) = epsq(i)
376 sminrl(i) = svm(i)
377 ENDDO
378 END IF ! unloading
379 END IF ! FUND > 0
380c-------------------------------------------------------------------------
381 uvar(1:nel,1) = emax(1:nel)
382 uvar(1:nel,2) = smax(1:nel)
383 uvar(1:nel,3) = eminrl(1:nel)
384 uvar(1:nel,4) = emaxrl(1:nel)
385 uvar(1:nel,5) = sminrl(1:nel)
386 uvar(1:nel,6) = smaxrl(1:nel)
387 uvar(1:nel,7) = zero
388 uvar(1:nel,10)= zero
389 DO ii=1,nindxc
390 i = indxc(ii)
391 uvar(i,10) = -one
392 END DO
393 DO i=1,nel
394 uvar(i,8) = sqrt((epsxx(i)**2 + epsyy(i)**2) / (one + nu21**2))
395 END DO
396c-----------
397 RETURN
398 END
interface_table_mod::table_vinterp
Definition table_mod.F:251
law119_membrane
subroutine law119_membrane(nel, nuparam, nuvar, uparam, uvar, gs, et, depsxx, depsyy, depsxy, epsxx, epsyy, epsxy, epsyz, epszx, sigoxx, sigoyy, sigoxy, sigoyz, sigozx, signxx, signyy, signxy, signyz, signzx, numtabl, itable, table)
Definition law119_membrane.F:41
max
#define max(a, b)
Definition macros.h:21
interface_table_mod
Definition table_mod.F:238
message_mod
Definition message_mod.F:1249
table_mod
Definition table_mod.F:112
table_mod::ttable
Definition table_mod.F:125