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,J,K,ITER,NORDER
77 my_real :: NU,RBULK,TENSCUT,GMAX,RVT,SUM,SUMDWDL,PARTP,EMAX,A11
78 my_real
79 . MU(5),AL(5),EVMA1(NEL,5),EVMA2(NEL,5),EVMA3(NEL,5),EVM(NEL,3),
80 . EIGV(NEL,3,2),TRAV(NEL),ROOTV(NEL),EVV(NEL,3),EV(NEL,3),
81 . RHO(NEL),DEZZ(NEL),DWDL(3),DDWDDL(3),RV(NEL),T(NEL,3),
82 . EPSZZ(NEL),SV(NEL,3),
83 . h30(100),h31(nel,100),h1(100),h10(100),
84 . h2(100),h20(100),h12(100),h120(100),
85 . cd1(nel),cd2(nel),cd12(nel),cd10(nel),
86 . cd20(nel),cd120(nel), dc3ev3(nel),c31(nel),c30(nel),kt3(nel)
87C=======================================================================
88 mu(1) = uparam(1)
89 mu(2) = uparam(2)
90 mu(3) = uparam(3)
91 mu(4) = uparam(4)
92 mu(5) = uparam(5)
93 al(1) = uparam(6)
94 al(2) = uparam(7)
95 al(3) = uparam(8)
96 al(4) = uparam(9)
97 al(5) = uparam(10)
98 rbulk = uparam(11)
99 tenscut= uparam(12)
100 nu = uparam(14)
101 norder = nint(uparam(18))
102C
103 gmax = zero
104 DO i=1,norder
105 gmax = gmax + mu(i)*al(i)
106 ENDDO
107C User variables initialisation
108 IF (time == zero .AND. isigi == 0) THEN
109 DO i=1,nel
110 DO j=1,nuvar
111 uvar(i,j) = zero
112 ENDDO
113 uvar(i,3) = one
114 ENDDO
115 ENDIF
116C principal stretch (def gradient eigenvalues)
117 DO i=1,nel
118 trav(i) = epsxx(i)+epsyy(i)
119 rootv(i) = sqrt((epsxx(i)-epsyy(i))*(epsxx(i)-epsyy(i))
120 . + epsxy(i)*epsxy(i))
121 evv(i,1) = half*(trav(i)+rootv(i))
122 evv(i,2) = half*(trav(i)-rootv(i))
123 evv(i,3) = zero
124 ENDDO
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
163C--------------------------------------
164C Newton method => Find EV(3) : T3(EV(3)) = 0
165C--------------------------------------
166 DO iter = 1,5
167! ----------------
168 DO i=1,nel
169 rv(i) = ev(i,1)*ev(i,2)*ev(i,3)
170c---- normalized stretch => unified compressible/uncompressible formution
171! RVT = RV(I)**(-THIRD)
172 IF(rv(i)> zero) THEN
173 rvt = exp((-third)*log(rv(i)))
174 ELSE
175 rvt = zero
176 ENDIF
177 evm(i,1) = ev(i,1)*rvt
178 evm(i,2) = ev(i,2)*rvt
179 evm(i,3) = ev(i,3)*rvt
180 ENDDO ! 1,NEL
181! ----------------
182C---- partial derivatives of strain energy
183! ----------------
184 DO j=1,5
185 DO i=1,nel
186 IF(evm(i,1)>zero) THEN
187 IF(al(j)/=zero) THEN
188 evma1(i,j) = mu(j) * exp(al(j)*log(evm(i,1)))
189 ELSE
190 evma1(i,j) = mu(j)
191 ENDIF
192 ELSE
193 evma1(i,j) = zero
194 ENDIF
195 IF(evm(i,2)>zero) THEN
196 IF(al(j)/=zero) THEN
197 evma2(i,j) = mu(j) * exp(al(j)*log(evm(i,2)))
198 ELSE
199 evma2(i,j) = mu(j)
200 ENDIF
201 ELSE
202 evma2(i,j) = zero
203 ENDIF
204 IF(evm(i,3)>zero) THEN
205 IF(al(j)/=zero) THEN
206 evma3(i,j) = mu(j) * exp(al(j)*log(evm(i,3)))
207 ELSE
208 evma3(i,j) = mu(j)
209 ENDIF
210 ELSE
211 evma3(i,j) = zero
212 ENDIF
213 ENDDO ! 1,NEL
214 ENDDO ! j=1,5
215! ----------------
216 DO i=1,nel
217 dwdl(1) = evma1(i,1)+evma1(i,2)+evma1(i,3)+evma1(i,4)+evma1(i,5)
218 dwdl(2) = evma2(i,1)+evma2(i,2)+evma2(i,3)+evma2(i,4)+evma2(i,5)
219 dwdl(3) = evma3(i,1)+evma3(i,2)+evma3(i,3)+evma3(i,4)+evma3(i,5)
220 sumdwdl = (dwdl(1)+dwdl(2)+dwdl(3))* third
221 partp = rbulk*(rv(i)- one)
222c------ ---
223c principal cauchy stress
224 t(i,1) = (dwdl(1) - sumdwdl) / rv(i) + partp
225 t(i,2) = (dwdl(2) - sumdwdl) / rv(i) + partp
226 t(i,3) = (dwdl(3) - sumdwdl) / rv(i) + partp
227c------ ---
228
229 kt3(i) = -third*(evma1(i,1)+evma1(i,2)+evma1(i,3)+evma1(i,4)+evma1(i,5))
230 . -third*(evma2(i,1)+evma2(i,2)+evma2(i,3)+evma2(i,4)+evma2(i,5))
231 . +two_third*(evma3(i,1)+evma3(i,2)+evma3(i,3)+evma3(i,4)+evma3(i,5))
232 kt3(i) =-ev(i,1)*ev(i,2)*kt3(i)/(rv(i)**2) + rbulk*ev(i,1)*ev(i,2)
233 kt3(i) = kt3(i)
234 . +(one_over_9*(al(1)*evma1(i,1)+al(2)*evma1(i,2)+al(3)*evma1(i,3)
235 . + al(4)*evma1(i,4)+al(5)*evma1(i,5)
236 . + al(1)*evma2(i,1)+al(2)*evma2(i,2) + al(3)*evma2(i,3)
237 . + al(4)*evma2(i,4)+al(5)*evma2(i,5)
238 . + four*(al(1)*evma3(i,1) + al(2)*evma3(i,2)
239 . + al(3)*evma3(i,3)
240 . + al(4)*evma3(i,4)+al(5)*evma3(i,5))))/ev(i,3)/rv(i)
241C
242 ev(i,3) = ev(i,3) - t(i,3)/kt3(i)
243 rv(i) = ev(i,1)*ev(i,2)*ev(i,3)
244 ENDDO ! 1,NEL
245 ENDDO ! ITER = 1,5
246C-----------------------
247 sv(1:nel,1) = zero
248 sv(1:nel,2) = zero
249 sv(1:nel,3) = zero
250C-----------------------
251 DO i=1,nel
252 uvar(i,3) = ev(i,3)
253 ENDDO
254C--------------------------------------
255c tension cut
256 DO i=1,nel
257 IF (off(i) /= zero .AND.
258 . (t(i,1) > abs(tenscut) .OR. t(i,2) > abs(tenscut))) THEN
259 t(i,1) = zero
260 t(i,2) = zero
261 t(i,3) = zero
262 off(i) = four_over_5
263 ENDIF
264 ENDDO
265C-------------------------------------------------------------
266C set sound speed & viscosity
267 ! DO I=1,NEL
268 ! DEZZ(I) =-NU/(ONE-NU)*(DEPSXX(I)+DEPSYY(I))
269 ! THKN(I) = THKN(I) + DEZZ(I)*THKLYL(I)
270 ! RHO(I) = RHO0(I)/RV(I)
271 ! SOUNDSP(I) = SQRT((TWO_THIRD*GMAX+RBULK)/RHO(I))
272 ! VISCMAX(I) = ZERO
273 ! ENDDO
274 IF (ismstr == 1 .OR. ismstr == 3 .OR. ismstr == 11) THEN ! engineering strain
275 DO i=1,nel
276 epszz(i) =ev(i,3) - one
277 uvar(i,3) = ev(i,3)
278 ENDDO
279 ELSEIF (ismstr == 10) THEN ! left gauchy-green strain
280 DO i=1,nel
281 epszz(i) =ev(i,3) - one
282 uvar(i,3) = ev(i,3)
283 ENDDO
284 ELSE ! true strain
285 DO i=1,nel
286 epszz(i) =log(ev(i,3))
287 uvar(i,3) = ev(i,3)
288 ENDDO
289 ENDIF
290 DO i=1,nel
291 rv(i) = ev(i,1)*ev(i,2)*ev(i,3)
292 dezz(i) =-nu/(one-nu)*(depsxx(i)+depsyy(i))
293 signxx(i) =eigv(i,1,1)*t(i,1)+eigv(i,1,2)*t(i,2) + sv(i,1)/rv(i)
294 signyy(i) =eigv(i,2,1)*t(i,1)+eigv(i,2,2)*t(i,2) + sv(i,2)/rv(i)
295 signxy(i) =eigv(i,3,1)*t(i,1)+eigv(i,3,2)*t(i,2) + sv(i,3)/rv(i)
296C
297 signyz(i) = sigoyz(i)+gs(i)*depsyz(i)
298 signzx(i) = sigozx(i)+gs(i)*depszx(i)
299 rho(i) = rho0(i)/rv(i)
300 thkn(i) = thkn(i) + dezz(i)*thklyl(i)*off(i)
301 viscmax(i)= zero
302C
303 emax = gmax*(one + nu)
304 a11 = emax/(one - nu**2)
305 soundsp(i)= sqrt(a11/rho0(i))
306 ENDDO
307C-----------
308 RETURN
309 END
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