OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps82c.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!|| sigeps82c ../engine/source/materials/mat/mat082/sigeps82c.F
25!||--- called by ------------------------------------------------------
26!|| mulawc ../engine/source/materials/mat_share/mulawc.F90
27!||--- calls -----------------------------------------------------
28!|| finter ../engine/source/tools/curve/finter.F
29!||====================================================================
30 SUBROUTINE sigeps82c(
31 1 NEL , NUPARAM, NUVAR , NFUNC , IFUNC ,
32 2 NPF , NPT0 , IPT , IFLAG ,
33 2 TF , TIME , TIMESTEP,
34 3 UPARAM , RHO0 ,AREA , EINT , THKLYL,
35 4 EPSPXX , EPSPYY , EPSPXY, EPSPYZ, EPSPZX,
36 5 DEPSXX , DEPSYY , DEPSXY, DEPSYZ, DEPSZX,
37 6 EPSXX , EPSYY , EPSXY , EPSYZ , EPSZX ,
38 7 SIGOXX , SIGOYY , SIGOXY, SIGOYZ, SIGOZX,
39 8 SIGNXX , SIGNYY , SIGNXY, SIGNYZ, SIGNZX,
40 9 SIGVXX , SIGVYY , SIGVXY, SIGVYZ, SIGVZX,
41 A SOUNDSP, VISCMAX, THKN , UVAR ,
42 B NGL , OFF , ISMSTR , IPM , GS )
43C-----------------------------------------------
44C I M P L I C I T T Y P E S
45C-----------------------------------------------
46#include "implicit_f.inc"
47C-----------------------------------------------
48C C O M M O N
49C-----------------------------------------------
50#include "param_c.inc"
51C----------------------------------------------------------------
52C I N P U T A R G U M E N T S
53C----------------------------------------------------------------
54 INTEGER NEL,NUPARAM, NUVAR,ISMSTR, IPM(NPROPMI,*),MAT(NEL) ,
55 . NPT0,IPT,IFLAG(*),NGL(NEL)
56 my_real
57 . TIME,TIMESTEP,UPARAM(NUPARAM),THKN(NEL),THKLYL(NEL),
58 . RHO0(NEL),AREA(NEL),EINT(NEL,2),GS(NEL),
59 . EPSPXX(NEL),EPSPYY(NEL),EPSPXY(NEL),EPSPYZ(NEL),EPSPZX(NEL),
60 . DEPSXX(NEL),DEPSYY(NEL),DEPSXY(NEL),DEPSYZ(NEL),DEPSZX(NEL),
61 . EPSXX (NEL),EPSYY (NEL),EPSXY (NEL),EPSYZ (NEL),EPSZX (NEL),
62 . SIGOXX(NEL),SIGOYY(NEL),SIGOXY(NEL),SIGOYZ(NEL),SIGOZX(NEL)
63C----------------------------------------------------------------
64C O U T P U T A R G U M E N T S
65C----------------------------------------------------------------
66 my_real
67 . signxx(nel),signyy(nel),signxy(nel),signyz(nel),signzx(nel),
68 . sigvxx(nel),sigvyy(nel),sigvxy(nel),sigvyz(nel),sigvzx(nel),
69 . soundsp(nel),viscmax(nel),et(nel)
70C----------------------------------------------------------------
71C I N P U T O U T P U T A R G U M E N T S
72C----------------------------------------------------------------
73 my_real
74 . uvar(nel,nuvar), off(nel)
75C----------------------------------------------------------------
76C VARIABLES FOR FUNCTION INTERPOLATION
77C----------------------------------------------------------------
78 INTEGER NPF(*), NFUNC, IFUNC(NFUNC)
79 my_real FINTER,FINTTE,TF(*),FINT2V
80 EXTERNAL FINTER,FINTTE
81C----------------------------------------------------------------
82C L O C A L V A R I B L E S
83C----------------------------------------------------------------
84 INTEGER I,K,ITER,NORDRE,I_K2
85 my_real
86 . nu,tenscut,gmax,rbulk,ffac,rvt,partt(nel),partp(nel),sum,dd,k2
87 my_real
88 . evv(nel,3),ev(nel,3),evm(nel,3),rv(nel),rho(nel),dwdl(nel,3),
89 . dezz(nel),eigv(nel,3,2),trav(nel),rootv(nel),t(nel,3)
90 my_real
91 . mu(100),al(100),d(100),evma1(nel,100),evma2(nel,100),evma3(nel,100)
92C=======================================================================
93C SET INITIAL MATERIAL CONSTANTS
94 gmax = zero
95 nordre = nint(uparam(1))
96 DO i= 1,nordre
97 mu(i) = uparam(1 + i)
98 al(i) = uparam(1 + nordre + i)
99 d(i) = uparam(1 + 2*nordre + i)
100 gmax = gmax + mu(i)
101 ENDDO
102 rbulk = two/d(1)
103 nu = (three*rbulk-two*gmax)/(two*gmax+six*rbulk)
104 IF (nu == half) nu = 0.495
105 IF (time == zero) uvar(1:nel,1)=one
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 rot matrix (eigenvectors)
116 DO i=1,nel
117 IF (abs(evv(i,2)-evv(i,1)) < em10) THEN
118 eigv(i,1,1) = one
119 eigv(i,2,1) = one
120 eigv(i,3,1) = zero
121 eigv(i,1,2) = zero
122 eigv(i,2,2) = zero
123 eigv(i,3,2) = zero
124 ELSE
125 eigv(i,1,1) = one/rootv(i)*(epsxx(i)-evv(i,2))
126 eigv(i,2,1) = one/rootv(i)*(epsyy(i)-evv(i,2))
127 eigv(i,1,2) = one/rootv(i)*(evv(i,1)-epsxx(i))
128 eigv(i,2,2) = one/rootv(i)*(evv(i,1)-epsyy(i))
129 eigv(i,3,1) = one/rootv(i)*(half*epsxy(i))
130 eigv(i,3,2) =-one/rootv(i)*(half*epsxy(i))
131 ENDIF
132 ENDDO
133C--- Strain definition
134 IF (ismstr == 1 .OR. ismstr == 3) THEN ! engineering strain
135 DO i=1,nel
136 ev(i,1) = evv(i,1) + one
137 ev(i,2) = evv(i,2) + one
138 ev(i,3) = uvar(i,1)
139 ENDDO
140 ELSE ! true strain
141 DO i=1,nel
142 ev(i,1) = exp(evv(i,1))
143 ev(i,2) = exp(evv(i,2))
144 ev(i,3) = uvar(i,1)
145 ENDDO
146 ENDIF
147C--------------------------------------
148C Newton method => Find root EV(3) : T3(EV(3)) = 0
149C--------------------------------------
150 DO iter = 1,3
151 rv(1:nel) = ev(1:nel,1)*ev(1:nel,2)*ev(1:nel,3)
152c---- normalized stretch => unified compressible/uncompressible formulation
153 DO i=1,nel
154 IF(rv(i)/=zero) THEN
155 ! RVT = RV(I)**(-THIRD)
156 rvt = rv(i)**(-third) ! EXP((-THIRD)*LOG(RV(I)))
157 ELSE
158 rvt = zero
159 ENDIF
160 evm(i,1) = ev(i,1)*rvt
161 evm(i,2) = ev(i,2)*rvt
162 evm(i,3) = ev(i,3)*rvt
163 ENDDO
164 DO k = 1,nordre
165 DO i=1,nel
166 IF(evm(i,1)/=zero) THEN
167 evma1(i,k) = evm(i,1)**al(k)
168 ELSE
169 evma1(i,k) = zero
170 ENDIF
171 IF(evm(i,2)/=zero) THEN
172 evma2(i,k) = evm(i,2)**al(k)
173 ELSE
174 evma2(i,k) = zero
175 ENDIF
176 IF(evm(i,3)/=zero) THEN
177 evma3(i,k) = evm(i,3)**al(k)
178 ELSE
179 evma3(i,k) = zero
180 ENDIF
181 ENDDO ! 1,NEL
182 ENDDO ! Nordre
183
184C---- partial derivatives of strain energy
185 partt(1:nel) =zero
186 DO k = 1,nordre
187 dd = two*mu(k)/al(k)
188 DO i=1,nel
189 sum = third*(evma1(i,k) + evma2(i,k) + evma3(i,k))
190 partt(i) = partt(i) + dd*(evma3(i,k) - sum)
191 ENDDO
192 ENDDO
193 partp(1:nel) =zero
194 DO k = 1,nordre
195 IF (d(k) /= zero) THEN
196 k2 = two*k
197 i_k2 = 2*k
198 dd = k2/d(k)
199 DO i=1,nel
200 partp(i) = partp(i) + dd*(rv(i)-one)**(i_k2-1)
201 ENDDO
202 ENDIF
203 ENDDO
204C--------------------------------
205 t(1:nel,3) = partt(1:nel)/rv(1:nel) + partp(1:nel)
206C--------------------------------
207c second derivative (d_sigma/d_lambda3)
208 partt(1:nel) =zero
209 DO k = 1,nordre
210 partt(1:nel) = partt(1:nel)
211 . + two*mu(k)*(evma1(1:nel,k)+evma2(1:nel,k)+four*evma3(1:nel,k))/nine
212 ENDDO
213 partp(1:nel) = zero
214 DO k = 1,nordre
215 IF (d(k) /= zero) THEN
216 k2 = two*k
217 i_k2 = 2*k
218 dd = k2*(k2-one)/d(k)
219! PARTP = PARTP + DD*(RV(I)-ONE)**(K2-TWO)
220 partp(1:nel) = partp(1:nel) + dd*(rv(1:nel)-one)**(i_k2-2)
221 ENDIF
222 ENDDO
223 partt(1:nel) = partt(1:nel) / rv(1:nel)+ partp(1:nel) - t(1:nel,3)
224 ev(1:nel,3) = ev(1:nel,3)*(one - t(1:nel,3)/partt(1:nel))
225 ENDDO ! ITER
226C----------------
227c Recalculate the principal stress
228 DO i=1,nel
229 rv(i) = ev(i,1)*ev(i,2)*ev(i,3)
230! RVT = RV(I)**(-THIRD)
231 IF(rv(i)/=zero) THEN
232 rvt = rv(i)**(-third) !EXP((-THIRD)*LOG(RV(I)))
233 ELSE
234 rvt = zero
235 ENDIF
236 evm(i,1) = ev(i,1)*rvt
237 evm(i,2) = ev(i,2)*rvt
238 evm(i,3) = ev(i,3)*rvt
239 ENDDO
240 DO k = 1,nordre
241 DO i=1,nel
242 IF(evm(i,1)/=zero) THEN
243 evma1(i,k) = evm(i,1)**al(k)
244 ELSE
245 evma1(i,k) = zero
246 ENDIF
247 IF(evm(i,2)/=zero) THEN
248 evma2(i,k) = evm(i,2)**al(k)
249 ELSE
250 evma2(i,k) = zero
251 ENDIF
252 IF(evm(i,3)>zero) THEN
253 evma3(i,k) = evm(i,3)**al(k)
254 ELSE
255 evma3(i,k) = zero
256 ENDIF
257 ENDDO ! 1:NEL
258 ENDDO ! Nordre
259 dwdl(1:nel,1)=zero
260 dwdl(1:nel,2)=zero
261 dwdl(1:nel,3)=zero
262 DO k = 1,nordre
263 dd = mu(k)/al(k)
264 DO i=1,nel
265 sum = third*(evma1(i,k) + evma2(i,k) + evma3(i,k))
266 dwdl(i,1) = dwdl(i,1) + dd*(evma1(i,k) - sum)
267 dwdl(i,2) = dwdl(i,2) + dd*(evma2(i,k) - sum)
268 dwdl(i,3) = dwdl(i,3) + dd*(evma3(i,k) - sum)
269 ENDDO
270 ENDDO
271 partp(1:nel) = zero
272 DO k = 1,nordre
273 IF (d(k) /= zero) THEN
274 k2 = two*k
275 i_k2 = 2*k
276 dd = k2/d(k)
277 partp(1:nel) = partp(1:nel) + dd*(rv(1:nel)-one)**(i_k2-1)
278 ENDIF
279 ENDDO
280C
281 t(1:nel,1) = two*dwdl(1:nel,1)/rv(1:nel) + partp(1:nel)
282 t(1:nel,2) = two*dwdl(1:nel,2)/rv(1:nel) + partp(1:nel)
283 t(1:nel,3) = two*dwdl(1:nel,3)/rv(1:nel) + partp(1:nel)
284C-------------------------------------------------------------
285 DO i=1,nel
286 uvar(i,1) = ev(i,3)
287 ENDDO
288C-------------------------------------------------------------
289C transform principal Cauchy stress to global directions
290 DO i=1,nel
291 signxx(i) = eigv(i,1,1)*t(i,1) + eigv(i,1,2)*t(i,2)
292 signyy(i) = eigv(i,2,1)*t(i,1) + eigv(i,2,2)*t(i,2)
293 signxy(i) = eigv(i,3,1)*t(i,1) + eigv(i,3,2)*t(i,2)
294 signyz(i) = sigoyz(i) + gs(i)*depsyz(i)
295 signzx(i) = sigozx(i) + gs(i)*depszx(i)
296 ENDDO
297C-------------------------------------------------------------
298C set thickness, sound speed & viscosity
299 DO i=1,nel
300 dezz(i) =-nu/(one-nu)*(depsxx(i)+depsyy(i))
301 thkn(i) = thkn(i) + dezz(i)*thklyl(i)
302 rho(i) = rho0(i)/rv(i)
303 soundsp(i) = sqrt((two_third*gmax+rbulk)/rho(i))
304 viscmax(i) = zero
305 ENDDO
306C-----------
307 RETURN
308 END
sigeps82c
subroutine sigeps82c(nel, nuparam, nuvar, nfunc, ifunc, npf, npt0, ipt, iflag, tf, 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, ipm, gs)
Definition sigeps82c.F:43