OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps69c.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!|| sigeps69c ../engine/source/materials/mat/mat069/sigeps69c.F
25!||--- called by ------------------------------------------------------
26!|| mulawc ../engine/source/materials/mat_share/mulawc.F90
27!||====================================================================
28 SUBROUTINE sigeps69c(
29 1 NEL , NUPARAM, NUVAR ,NPT0 , ILAYER ,
30 2 TIME , TIMESTEP, UPARAM, RHO0 ,
31 3 AREA , EINT , THKLYL,
32 4 EPSPXX , EPSPYY , EPSPXY, EPSPYZ, EPSPZX,
33 5 DEPSXX , DEPSYY , DEPSXY, DEPSYZ, DEPSZX,
34 6 EPSXX , EPSYY , EPSXY , EPSYZ , EPSZX ,
35 7 SIGOXX , SIGOYY , SIGOXY, SIGOYZ, SIGOZX,
36 8 SIGNXX , SIGNYY , SIGNXY, SIGNYZ, SIGNZX,
37 9 SIGVXX , SIGVYY , SIGVXY, SIGVYZ, SIGVZX,
38 A SOUNDSP, VISCMAX, THKN , UVAR , NGL ,
39 B OFF , ISMSTR , GS )
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 "param_c.inc"
48#include "com01_c.inc"
49C----------------------------------------------------------------
50C I N P U T A R G U M E N T S
51C----------------------------------------------------------------
52 INTEGER NEL,NUPARAM,NUVAR,ISMSTR,NPT0,ILAYER
53 INTEGER NGL(NEL)
54 my_real :: TIME,TIMESTEP
55 my_real
56 . UPARAM(*),THKN(NEL),THKLYL(NEL),
57 . RHO0(NEL),AREA(NEL),EINT(NEL,2),GS(NEL),
58 . EPSPXX(NEL),EPSPYY(NEL),EPSPXY(NEL),EPSPYZ(NEL),EPSPZX(NEL),
59 . DEPSXX(NEL),DEPSYY(NEL),DEPSXY(NEL),DEPSYZ(NEL),DEPSZX(NEL),
60 . EPSXX (NEL),EPSYY (NEL),EPSXY (NEL),EPSYZ (NEL),EPSZX (NEL),
61 . SIGOXX(NEL),SIGOYY(NEL),SIGOXY(NEL),SIGOYZ(NEL),SIGOZX(NEL)
62C----------------------------------------------------------------
63C O U T P U T A R G U M E N T S
64C----------------------------------------------------------------
65 my_real
66 . signxx(nel),signyy(nel),signxy(nel),signyz(nel),signzx(nel),
67 . sigvxx(nel),sigvyy(nel),sigvxy(nel),sigvyz(nel),sigvzx(nel),
68 . soundsp(nel),viscmax(nel)
69C----------------------------------------------------------------
70C I N P U T O U T P U T A R G U M E N T S
71C----------------------------------------------------------------
72 my_real :: uvar(nel,nuvar), off(nel)
73C----------------------------------------------------------------
74C L O C A L V A R I B L E S
75C----------------------------------------------------------------
76 INTEGER :: I,ITER,NORDER,NITER,II
77 my_real :: NU,RBULK,TENSCUT,GMAX,SUM,SUMDWDL,PARTP,EMAX,A11
78 my_real :: AMAX
79 my_real ,DIMENSION(3) :: LAM_AL
80 my_real ,DIMENSION(5) :: MU,AL
81 my_real ,DIMENSION(NEL) :: RVT,GTMAX,DLAM3,KT3,RHO,KIR3
82 my_real ,DIMENSION(NEL) :: INVRV,DEZZ,RV,TRAV,ROOTV
83 my_real ,DIMENSION(NEL,3) :: t,evv,ev,evm,cii,s_ldwdl
84 my_real ,DIMENSION(NEL,3,2) :: eigv(nel,3,2)
85C=======================================================================
86 mu(1) = uparam(1)
87 mu(2) = uparam(2)
88 mu(3) = uparam(3)
89 mu(4) = uparam(4)
90 mu(5) = uparam(5)
91 al(1) = uparam(6)
92 al(2) = uparam(7)
93 al(3) = uparam(8)
94 al(4) = uparam(9)
95 al(5) = uparam(10)
96 rbulk = uparam(11)
97 tenscut= uparam(12)
98 nu = uparam(14)
99 norder = nint(uparam(18))
100C
101 gmax = zero
102 DO i=1,norder
103 gmax = gmax + mu(i)*al(i)
104 ENDDO
105! initialization of uvar(,3) in Starter
106C principal stretch (def gradient eigenvalues)
107 DO i=1,nel
108 trav(i) = epsxx(i)+epsyy(i)
109 rootv(i) = sqrt((epsxx(i)-epsyy(i))*(epsxx(i)-epsyy(i))
110 . + epsxy(i)*epsxy(i))
111 evv(i,1) = half*(trav(i)+rootv(i))
112 evv(i,2) = half*(trav(i)-rootv(i))
113 evv(i,3) = zero
114 ENDDO
115C-- avoid NaN---------
116 IF (ismstr == 10) THEN
117 DO i=1,nel
118 IF (min(evv(i,1),evv(i,2)) <= -one) THEN
119 evv(i,1) = zero
120 evv(i,2) = zero
121 off(i) = four_over_5
122 END IF
123 ENDDO
124 END IF
125C rot matrix (eigenvectors)
126 DO i=1,nel
127 IF(abs(evv(i,2)-evv(i,1))<em10) THEN
128 eigv(i,1,1) = one
129 eigv(i,2,1) = one
130 eigv(i,3,1) = zero
131 eigv(i,1,2) = zero
132 eigv(i,2,2) = zero
133 eigv(i,3,2) = zero
134 ELSE
135 eigv(i,1,1) = (epsxx(i)-evv(i,2)) /rootv(i)
136 eigv(i,2,1) = (epsyy(i)-evv(i,2)) /rootv(i)
137 eigv(i,1,2) = (evv(i,1)-epsxx(i)) /rootv(i)
138 eigv(i,2,2) = (evv(i,1)-epsyy(i)) /rootv(i)
139 eigv(i,3,1) = (half*epsxy(i)) /rootv(i)
140 eigv(i,3,2) =-(half*epsxy(i)) /rootv(i)
141 ENDIF
142 ENDDO
143C Strain definition
144 IF (ismstr == 1 .OR. ismstr == 3 .OR. ismstr == 11) THEN ! engineering strain
145 DO i=1,nel
146 ev(i,1)=evv(i,1)+ one
147 ev(i,2)=evv(i,2)+ one
148 ev(i,3)=uvar(i,3)
149 ENDDO
150 ELSEIF(ismstr == 10) THEN
151 DO i=1,nel
152 ev(i,1)=sqrt(evv(i,1)+ one)
153 ev(i,2)=sqrt(evv(i,2)+ one)
154 ev(i,3)=one/ev(i,1)/ev(i,2)
155 ENDDO
156 ELSE ! true strain
157 DO i=1,nel
158 ev(i,1)=exp(evv(i,1))
159 ev(i,2)=exp(evv(i,2))
160 ev(i,3)=uvar(i,3)
161 ENDDO
162 ENDIF
163 niter = 4
164!--------------------------------------
165! Newton method => Find EV(3) : Kirchoff J*T3(EV(3)) = 0
166!--------------------------------------
167 dlam3(1:nel) =zero
168 DO iter = 1,niter
169 ev(1:nel,3) = uvar(1:nel,3)+dlam3(1:nel) ! initial value takes lamda3(t)
170 DO i=1,nel
171 rv(i) = ev(i,1)*ev(i,2)*ev(i,3)
172 rvt(i) = exp( (-third)* log(rv(i)) )
173 evm(i,1:3) = ev(i,1:3)*rvt(i)
174 ENDDO ! 1,NEL
175 kir3(1:nel) = zero
176 kt3(1:nel) = zero
177 DO ii = 1,norder
178 IF (mu(ii)*al(ii) /= zero) THEN
179 DO i=1,nel
180 IF (off(i)==zero.OR.off(i)==four_over_5) cycle
181 lam_al(1:3) = exp(al(ii)*log(evm(i,1:3)))
182 sumdwdl = third*(lam_al(1)+lam_al(2)+lam_al(3))
183 sum = mu(ii)*(lam_al(3)-sumdwdl)
184 kir3(i) = kir3(i) + sum
185 kt3(i) = kt3(i) + al(ii)*sum
186 ENDDO
187 ENDIF
188 ENDDO
189 DO i=1,nel
190 IF (off(i)==zero.OR.off(i)==four_over_5) cycle
191 partp = rbulk*(rv(i)- one)
192 t(i,3)= kir3(i) + partp*rv(i) !Kirchoff
193 kt3(i)= two_third*kt3(i)/ev(i,3)+ rbulk*(two*rv(i)-one)*ev(i,1)*ev(i,2)
194 IF (kt3(i)>em20) dlam3(i) = dlam3(i) -t(i,3)/kt3(i)
195 ENDDO
196 END DO ! ITER = 1,NITER
197 ev(1:nel,3) = uvar(1:nel,3)+dlam3(1:nel)
198 dezz(1:nel) = log(one+dlam3(1:nel)/uvar(1:nel,3))
199 uvar(1:nel,3) = ev(1:nel,3)
200! compute T1,T2 Cauchy stress
201 DO i=1,nel
202 rv(i) = ev(i,1)*ev(i,2)*ev(i,3)
203 rvt(i) = exp((-third)*log(rv(i))) ! -> J^(-1/3)
204 evm(i,1:3) = ev(i,1:3)*rvt(i)
205 invrv(i) = one / rv(i)
206 END DO
207 s_ldwdl(1:nel,1:3) = zero
208 DO ii = 1,norder
209 IF (mu(ii)*al(ii) /= zero) THEN
210 DO i=1,nel
211 lam_al(1:3) = exp(al(ii)*log(evm(i,1:3)))
212 s_ldwdl(i,1:3) = s_ldwdl(i,1:3) + mu(ii)*lam_al(1:3)
213 ENDDO
214 ENDIF
215 ENDDO
216 DO i=1,nel
217 sumdwdl = (s_ldwdl(i,1) + s_ldwdl(i,2) + s_ldwdl(i,3)) * third
218 partp = rbulk*(rv(i)- one)
219 t(i,1) = (s_ldwdl(i,1) - sumdwdl) *invrv(i) + partp
220 t(i,2) = (s_ldwdl(i,2) - sumdwdl) *invrv(i) + partp
221 ENDDO
222C--------------------------------------
223c tension cut
224 DO i=1,nel
225 IF (off(i) /= zero .AND.
226 . (t(i,1) > abs(tenscut) .OR. t(i,2) > abs(tenscut))) THEN
227 t(i,1) = zero
228 t(i,2) = zero
229 t(i,3) = zero
230 off(i) = four_over_5
231 ENDIF
232 ENDDO
233! new gt
234 cii(1:nel,1:3) = zero
235 DO ii = 1,norder
236 IF (mu(ii)*al(ii) /= zero) THEN
237 DO i=1,nel
238 lam_al(1:3) = exp(al(ii)*log(evm(i,1:3)))
239 amax = third*(lam_al(1)+lam_al(2)+lam_al(3))
240 cii(i,1:3) = cii(i,1:3) + mu(ii)*al(ii)*(lam_al(1:3)+amax)
241 ENDDO
242 ENDIF
243 ENDDO
244 DO i = 1,nel
245 gtmax(i) = half*max(cii(i,1),cii(i,2),cii(i,3))
246 ENDDO
247C-------------------------------------------------------------
248 DO i=1,nel
249 signxx(i) = eigv(i,1,1)*t(i,1) + eigv(i,1,2)*t(i,2)
250 signyy(i) = eigv(i,2,1)*t(i,1) + eigv(i,2,2)*t(i,2)
251 signxy(i) = eigv(i,3,1)*t(i,1) + eigv(i,3,2)*t(i,2)
252!
253 signyz(i) = sigoyz(i)+gs(i)*depsyz(i)
254 signzx(i) = sigozx(i)+gs(i)*depszx(i)
255 rho(i) = rho0(i)*invrv(i)
256 thkn(i) = thkn(i) + dezz(i)*thklyl(i)*off(i)
257 viscmax(i)= zero
258!
259 emax = max(gmax,gtmax(i))*(one + nu)
260 a11 = emax/(one - nu**2)
261 soundsp(i)= sqrt(a11/rho0(i))
262 ENDDO
263C-----------
264 RETURN
265 END
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
sigeps69c
subroutine sigeps69c(nel, nuparam, nuvar, npt0, ilayer, time, timestep, uparam, rho0, area, eint, thklyl, epspxx, epspyy, epspxy, epspyz, epspzx, depsxx, depsyy, depsxy, depsyz, depszx, epsxx, epsyy, epsxy, epsyz, epszx, sigoxx, sigoyy, sigoxy, sigoyz, sigozx, signxx, signyy, signxy, signyz, signzx, sigvxx, sigvyy, sigvxy, sigvyz, sigvzx, soundsp, viscmax, thkn, uvar, ngl, off, ismstr, gs)
Definition sigeps69c.F:40