OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps42c.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!|| sigeps42c ../engine/source/materials/mat/mat042/sigeps42c.F
25!||--- called by ------------------------------------------------------
26!|| mulawc ../engine/source/materials/mat_share/mulawc.F90
27!||====================================================================
28 SUBROUTINE sigeps42c(
29 1 NEL , NUPARAM ,NIPARAM, NUVAR , ISMSTR ,
30 2 TIME , TIMESTEP,UPARAM , IPARAM , RHO0 ,
31 3 DEPSXX , DEPSYY , DEPSXY , DEPSYZ , DEPSZX ,
32 4 EPSXX , EPSYY , EPSXY , THKN , THKLYL ,
33 5 SIGNXX , SIGNYY , SIGNXY , SIGNYZ , SIGNZX ,
34 6 SIGOYZ , SIGOZX , SOUNDSP, GS , UVAR ,
35 7 OFF )
36C-----------------------------------------------
37C I M P L I C I T T Y P E S
38C-----------------------------------------------
39#include "implicit_f.inc"
40C-----------------------------------------------
41C C O M M O N
42C-----------------------------------------------
43#include "param_c.inc"
44#include "com01_c.inc"
45C----------------------------------------------------------------
46C I N P U T A R G U M E N T S
47C----------------------------------------------------------------
48 INTEGER :: NEL
49 INTEGER :: NUPARAM
50 INTEGER :: NIPARAM
51 INTEGER :: NUVAR
52 INTEGER :: ISMSTR
53 INTEGER :: IPARAM(NIPARAM)
54 my_real :: UPARAM(NUPARAM)
55 my_real :: time,timestep
56 my_real ,DIMENSION(NEL) :: thkn,thklyl,rho0,gs,
57 . depsxx,depsyy,depsxy,depsyz,depszx,epsxx,epsyy,epsxy
58 my_real ,DIMENSION(NEL) ,INTENT(IN) :: sigoyz,sigozx
59C----------------------------------------------------------------
60C O U T P U T A R G U M E N T S
61C----------------------------------------------------------------
62 my_real ,DIMENSION(NEL) :: signxx,signyy,signxy,signyz,signzx,soundsp
63C----------------------------------------------------------------
64C I N P U T O U T P U T A R G U M E N T S
65C----------------------------------------------------------------
66 my_real :: uvar(nel,nuvar), off(nel)
67C----------------------------------------------------------------
68C L O C A L V A R I B L E S
69C----------------------------------------------------------------
70 INTEGER :: I,II,NPRONY,NORDER,ITER,JNV,NITER
71 my_real :: SUM,FAC,FSCAL,TENSCUT,SUMDWDL,PARTP,AMAX
72 my_real :: EMAX,GMAX,GVMAX,NU,RBULK,A11
73 my_real ,DIMENSION(3) :: lam_al
74 my_real ,DIMENSION(10) :: mu,al
75 my_real ,DIMENSION(NEL) :: rvt,gtmax,dlam3
76 my_real ,DIMENSION(NEL) :: invrv,dezz,rv,trav,rootv
77 my_real ,DIMENSION(NEL) :: kt3,rho,kir3
78 my_real ,DIMENSION(NEL,3) :: t,sv,evv,ev,evm,cii,s_ldwdl
79 my_real ,DIMENSION(NEL,3,2) :: eigv(nel,3,2)
80 my_real :: dav,h0(4),depszz,taux
81 my_real ,DIMENSION(NEL) :: a1,a2
82 my_real ,DIMENSION(NEL,4) :: devs
83 my_real, DIMENSION(:) , ALLOCATABLE :: gi,beta
84 my_real, DIMENSION(:,:) , ALLOCATABLE :: aa,bb,h11
85C=======================================================================
86C SET INITIAL MATERIAL CONSTANTS
87
88 norder = iparam(1)
89 nprony = iparam(2)
90 IF (nprony>0) THEN
91 ALLOCATE(aa(nel,nprony))
92 ALLOCATE(bb(nel,nprony))
93 ALLOCATE(h11(6,nprony))
94 ALLOCATE(gi(nprony))
95 ALLOCATE(beta(nprony))
96 ENDIF
97 sv(1:nel,1:3) = zero
98!
99 DO i=1,norder
100 mu(i) = uparam(i)
101 al(i) = uparam(i+10)
102 ENDDO
103 rbulk = uparam(21)
104 nu = uparam(22)
105 tenscut= uparam(23)
106 fscal = uparam(24)
107! ------------------
108 gmax = zero
109 gvmax = zero
110 DO i=1,nprony
111 gi(i) = uparam(24 + i)
112 taux = uparam(24 + nprony + i)
113 gvmax = gvmax + gi(i)
114 beta(i) = one/taux
115 ENDDO
116 DO i=1,norder
117 gmax = gmax + mu(i)*al(i)
118 ENDDO
119 gmax = gmax + gvmax
120C
121 niter = 4
122C principal stretch (def gradient eigenvalues)
123 DO i=1,nel
124 trav(i) = epsxx(i)+epsyy(i)
125 rootv(i) = sqrt((epsxx(i)-epsyy(i))*(epsxx(i)-epsyy(i))
126 . + epsxy(i)*epsxy(i))
127 evv(i,1) = half*(trav(i)+rootv(i))
128 evv(i,2) = half*(trav(i)-rootv(i))
129 evv(i,3) = zero
130 ENDDO
131C-- avoid NaN---------
132 IF (ismstr == 10) THEN
133 DO i=1,nel
134 IF (min(evv(i,1),evv(i,2)) <= -one) THEN
135 evv(i,1) = zero
136 evv(i,2) = zero
137 off(i) = four_over_5
138 END IF
139 ENDDO
140 END IF
141C rot matrix (eigenvectors)
142 DO i=1,nel
143 IF (abs(evv(i,2)-evv(i,1)) < em10) THEN
144 eigv(i,1,1)=one
145 eigv(i,2,1)=one
146 eigv(i,3,1)=zero
147 eigv(i,1,2)=zero
148 eigv(i,2,2)=zero
149 eigv(i,3,2)=zero
150 ELSE
151 invrv(i) = one / rootv(i)
152 eigv(i,1,1) = (epsxx(i)-evv(i,2)) * invrv(i)
153 eigv(i,2,1) = (epsyy(i)-evv(i,2)) * invrv(i)
154 eigv(i,3,1) = (half*epsxy(i)) * invrv(i)
155 eigv(i,1,2) = (evv(i,1)-epsxx(i)) * invrv(i)
156 eigv(i,2,2) = (evv(i,1)-epsyy(i)) * invrv(i)
157 eigv(i,3,2) =-(half*epsxy(i)) * invrv(i)
158 ENDIF
159 ENDDO
160C Strain definition
161 IF (ismstr == 1 .OR. ismstr == 3 .OR. ismstr == 11) THEN ! engineering strain
162 DO i=1,nel
163 ev(i,1)=evv(i,1)+ one
164 ev(i,2)=evv(i,2)+ one
165 ev(i,3)=one/ev(i,1)/ev(i,2)
166 ENDDO
167 ELSEIF(ismstr == 10) THEN
168 DO i=1,nel
169 ev(i,1)=sqrt(evv(i,1)+ one)
170 ev(i,2)=sqrt(evv(i,2)+ one)
171 ev(i,3)=one/ev(i,1)/ev(i,2)
172 ENDDO
173 ELSE ! true strain
174 DO i=1,nel
175 ev(i,1)=exp(evv(i,1))
176 ev(i,2)=exp(evv(i,2))
177 ev(i,3)=one/ev(i,1)/ev(i,2)
178 ENDDO
179 ENDIF
180 DO i=1,nel
181 IF (off(i)==zero.OR.off(i)==four_over_5) ev(i,1:3)=one
182 ENDDO
183! old iterartive Pnony removed, new formulation uncoupling w/ Prony
184!--------------------------------------
185! Newton method => Find EV(3) : Kirchoff J*T3(EV(3)) = 0
186!--------------------------------------
187 dlam3(1:nel) =zero
188 DO iter = 1,niter
189 ev(1:nel,3) = uvar(1:nel,3)+dlam3(1:nel) ! initial value takes lamda3(t)
190 DO i=1,nel
191 rv(i) = ev(i,1)*ev(i,2)*ev(i,3)
192 rvt(i) = exp( (-third)* log(rv(i)) )
193 evm(i,1:3) = ev(i,1:3)*rvt(i)
194 ENDDO ! 1,NEL
195 kir3(1:nel) = zero
196 kt3(1:nel) = zero
197 DO ii = 1,norder
198 IF (mu(ii)*al(ii) /= zero) THEN
199 DO i=1,nel
200 IF (off(i)==zero.OR.off(i)==four_over_5) cycle
201 lam_al(1:3) = exp(al(ii)*log(evm(i,1:3)))
202 sumdwdl = third*(lam_al(1)+lam_al(2)+lam_al(3))
203 sum = mu(ii)*(lam_al(3)-sumdwdl)
204 kir3(i) = kir3(i) + sum
205 kt3(i) = kt3(i) + al(ii)*sum
206 ENDDO
207 ENDIF
208 ENDDO
209 DO i=1,nel
210 IF (off(i)==zero.OR.off(i)==four_over_5) cycle
211 partp = rbulk*(rv(i)- one)
212 t(i,3)= kir3(i) + partp*rv(i) !Kirchoff
213 kt3(i)= two_third*kt3(i)/ev(i,3)+ rbulk*(two*rv(i)-one)*ev(i,1)*ev(i,2)
214 IF (kt3(i)>em20) dlam3(i) = dlam3(i) -t(i,3)/kt3(i)
215 ENDDO
216 END DO ! ITER = 1,NITER
217 ev(1:nel,3) = uvar(1:nel,3)+dlam3(1:nel) ! initial value takes lamda3(t): looking for (t+dt)
218 dezz(1:nel) = log(one+dlam3(1:nel)/uvar(1:nel,3))
219 uvar(1:nel,3) = ev(1:nel,3)
220! compute t1,t2 cauchy stress
221 DO i=1,nel
222 rv(i) = ev(i,1)*ev(i,2)*ev(i,3)
223 rvt(i) = exp((-third)*log(rv(i))) ! -> J^(-1/3)
224 evm(i,1:3) = ev(i,1:3)*rvt(i)
225 invrv(i) = one / rv(i)
226 END DO
227 s_ldwdl(1:nel,1:3) = zero
228 DO ii = 1,norder
229 IF (mu(ii)*al(ii) /= zero) THEN
230 DO i=1,nel
231 lam_al(1:3) = exp(al(ii)*log(evm(i,1:3)))
232 s_ldwdl(i,1:3) = s_ldwdl(i,1:3) + mu(ii)*lam_al(1:3)
233 ENDDO
234 ENDIF
235 ENDDO
236 DO i=1,nel
237 sumdwdl = (s_ldwdl(i,1) + s_ldwdl(i,2) + s_ldwdl(i,3)) * third
238 partp = rbulk*(rv(i)- one)
239 t(i,1) = (s_ldwdl(i,1) - sumdwdl) *invrv(i) + partp
240 t(i,2) = (s_ldwdl(i,2) - sumdwdl) *invrv(i) + partp
241 ENDDO
242 IF (nprony > 0 ) THEN ! Prony same as /VISC/PRONY
243 a1(1:nel) = zero
244 a2(1:nel) = zero
245 jnv = 12
246 DO ii=1,nprony
247 DO i=1,nel
248 aa(i,ii) = exp(-beta(ii)*timestep)
249 bb(i,ii) = gi(ii)*exp(-half*beta(ii)*timestep) !bb/dt, and use directly despsij
250! for computing DEPSZZ from (sigzz=0)
251 h0(3) = uvar(i,jnv + (ii - 1)*4 + 3)
252 a1(i) = a1(i) + aa(i,ii)*h0(3)
253 a2(i) = a2(i) + bb(i,ii)
254 END DO
255 END DO
256 DO i=1,nel
257 fac = one/max(em20,two_third*a2(i))
258 depszz = -a1(i)*fac + half*(depsxx(i) + depsyy(i))
259 dav = third*(depsxx(i) + depsyy(i) + depszz)
260 devs(i,1) = depsxx(i) - dav
261 devs(i,2) = depsyy(i) - dav
262 devs(i,3) = depszz - dav
263 devs(i,4) = half*depsxy(i)
264 dezz(i) = dezz(i) + depszz
265 END DO
266!
267 DO ii= 1,nprony
268 DO i=1,nel
269 h0(1:4) = uvar(i,jnv + (ii - 1)*4 +1:4)
270 h11(1:4,ii) = aa(i,ii)*h0(1:4) + bb(i,ii)*devs(i,1:4)
271 sv(i,1:2) = sv(i,1:2) + h11(1:2,ii)
272 sv(i,3) = sv(i,3) + h11(4,ii)
273 uvar(i,jnv + (ii - 1)*4 +1:4) = h11(1:4,ii)
274 END DO
275 END DO
276 END IF ! PRONY
277!-------------------------------------------------------------
278! tension cut
279 DO i=1,nel
280 IF (off(i) /= zero .AND.
281 . (t(i,1) > abs(tenscut) .OR. t(i,2) > abs(tenscut))) THEN
282 t(i,1) = zero
283 t(i,2) = zero
284 t(i,3) = zero
285 off(i) = four_over_5
286 ENDIF
287 ENDDO
288! new gt
289 cii(1:nel,1:3) = zero
290 DO ii = 1,norder
291 IF (mu(ii)*al(ii) /= zero) THEN
292 DO i=1,nel
293 lam_al(1:3) = exp(al(ii)*log(evm(i,1:3)))
294 amax = third*(lam_al(1)+lam_al(2)+lam_al(3))
295 cii(i,1:3) = cii(i,1:3) + mu(ii)*al(ii)*(lam_al(1:3)+amax)
296 ENDDO
297 ENDIF
298 ENDDO
299 DO i = 1,nel
300 gtmax(i) = gvmax+half*max(cii(i,1),cii(i,2),cii(i,3))
301 ENDDO
302!
303 DO i=1,nel
304 signxx(i) = eigv(i,1,1)*t(i,1) + eigv(i,1,2)*t(i,2) + sv(i,1)
305 signyy(i) = eigv(i,2,1)*t(i,1) + eigv(i,2,2)*t(i,2) + sv(i,2)
306 signxy(i) = eigv(i,3,1)*t(i,1) + eigv(i,3,2)*t(i,2) + sv(i,3)
307! transverse shear can't use visco(Prony) as /VISC due to incremental formulation
308 signyz(i) = sigoyz(i)+gs(i)*depsyz(i)
309 signzx(i) = sigozx(i)+gs(i)*depszx(i)
310 rho(i) = rho0(i)*invrv(i)
311 thkn(i) = thkn(i) + dezz(i)*thklyl(i)*off(i)
312
313!
314 emax = max(gmax,gtmax(i))*(one + nu)
315 a11 = emax/(one - nu**2)
316 soundsp(i)= sqrt(a11/rho0(i))
317 ENDDO
318 IF (nprony>0) THEN
319 DEALLOCATE(aa)
320 DEALLOCATE(bb)
321 DEALLOCATE(h11)
322 DEALLOCATE(gi)
323 DEALLOCATE(beta)
324 ENDIF
325C-----------
326 RETURN
327 END
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
sigeps42c
subroutine sigeps42c(nel, nuparam, niparam, nuvar, ismstr, time, timestep, uparam, iparam, rho0, depsxx, depsyy, depsxy, depsyz, depszx, epsxx, epsyy, epsxy, thkn, thklyl, signxx, signyy, signxy, signyz, signzx, sigoyz, sigozx, soundsp, gs, uvar, off)
Definition sigeps42c.F:36