OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps42.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!|| sigeps42 ../starter/source/materials/mat/mat042/sigeps42.F
25!||--- called by ------------------------------------------------------
26!|| mulaw ../starter/source/materials/mat_share/mulaw.F
27!||--- calls -----------------------------------------------------
28!|| finter ../starter/source/tools/curve/finter.F
29!|| valpvec ../starter/source/materials/tools/matrix.F
30!|| valpvecdp ../starter/source/materials/tools/matrix.F
31!||--- uses -----------------------------------------------------
32!||====================================================================
33 SUBROUTINE sigeps42(
34 1 MAT_PARAM,NEL ,NUVAR ,NFUNC ,IFUNC ,NPF ,
35 2 TF ,TIME ,TIMESTEP ,RHO0 ,RHO ,
36 3 EPSPXX ,EPSPYY ,EPSPZZ ,EPSPXY ,EPSPYZ ,EPSPZX ,
37 4 DEPSXX ,DEPSYY ,DEPSZZ ,DEPSXY ,DEPSYZ ,DEPSZX ,
38 5 EPSXX ,EPSYY ,EPSZZ ,EPSXY ,EPSYZ ,EPSZX ,
39 6 SIGOXX ,SIGOYY ,SIGOZZ ,SIGOXY ,SIGOYZ ,SIGOZX ,
40 7 SIGNXX ,SIGNYY ,SIGNZZ ,SIGNXY ,SIGNYZ ,SIGNZX ,
41 8 SOUNDSP ,VISCMAX ,UVAR ,OFF ,WXXDT ,WYYDT ,
42 A WZZDT ,ISMSTR ,MFXX ,MFXY ,MFXZ ,MFYX ,
43 B MFYY ,MFYZ ,MFZX ,MFZY ,MFZZ )
44C-----------------------------------------------
45C M o d u l e s
46C-----------------------------------------------
47 USE matparam_def_mod
48C-----------------------------------------------
49C I M P L I C I T T Y P E S
50C-----------------------------------------------
51#include "implicit_f.inc"
52C-----------------------------------------------
53C G L O B A L P A R A M E T E R S
54C-----------------------------------------------
55#include "mvsiz_p.inc"
56C-----------------------------------------------
57C C O M M O N
58C-----------------------------------------------
59#include "scr05_c.inc"
60#include "tabsiz_c.inc"
61C----------------------------------------------------------------
62C I N P U T A R G U M E N T S
63C----------------------------------------------------------------
64 INTEGER, INTENT(IN) :: NEL,NUVAR,ISMSTR
65 my_real, INTENT(IN) ::
66 . TIME , TIMESTEP ,
67 . RHO (NEL), RHO0 (NEL),
68 . EPSPXX(NEL), EPSPYY(NEL), EPSPZZ(NEL),
69 . EPSPXY(NEL), EPSPYZ(NEL), EPSPZX(NEL),
70 . DEPSXX(NEL), DEPSYY(NEL), DEPSZZ(NEL),
71 . DEPSXY(NEL), DEPSYZ(NEL), DEPSZX(NEL),
72 . EPSXX (NEL), EPSYY (NEL), EPSZZ (NEL),
73 . EPSXY (NEL), EPSYZ (NEL), EPSZX (NEL),
74 . sigoxx(nel), sigoyy(nel), sigozz(nel),
75 . sigoxy(nel), sigoyz(nel), sigozx(nel),
76 . mfxx(nel) , mfxy(nel), mfxz(nel),
77 . mfyx(nel) , mfyy(nel), mfyz(nel),
78 . mfzx(nel) , mfzy(nel), mfzz(nel),
79 . wzzdt(mvsiz),wyydt(mvsiz),wxxdt(mvsiz)
80 TYPE(matparam_struct_) ,INTENT(INOUT) :: MAT_PARAM
81C----------------------------------------------------------------
82C O U T P U T A R G U M E N T S
83C----------------------------------------------------------------
84 my_real, INTENT(OUT) ::
85 . SIGNXX (NEL), SIGNYY (NEL), SIGNZZ(NEL),
86 . SIGNXY (NEL), SIGNYZ (NEL), SIGNZX(NEL),
87 . SOUNDSP(NEL), VISCMAX(NEL)
88C----------------------------------------------------------------
89C I N P U T O U T P U T A R G U M E N T S
90C----------------------------------------------------------------
91 my_real, INTENT(INOUT) ::
92 . uvar(nel,nuvar), off(nel)
93C----------------------------------------------------------------
94C VARIABLES FOR FUNCTION INTERPOLATION
95C----------------------------------------------------------------
96 INTEGER NPF(SNPC),NFUNC,IFUNC(NFUNC)
97 my_real FINTER,TF(STF)
98 EXTERNAL FINTER
99C----------------------------------------------------------------
100C L O C A L V A R I B L E S
101C----------------------------------------------------------------
102 INTEGER I,J,K,KFP,IFORM,NORDER
103 my_real :: TENSIONCUT,GMAX,FFAC
104 my_real, DIMENSION(10) :: MU,AL
105 my_real
106 . sigprv(3,mvsiz),ev(3,mvsiz),evm(3,mvsiz),dwdl(3,mvsiz),
107 . t1(3,mvsiz),sumdwdl(mvsiz),rv(mvsiz),p(mvsiz),
108 . dwdrv(mvsiz),rbulk,dpdmu,av(6,mvsiz),evv(3,mvsiz),
109 . dirprv(3,3,mvsiz),a(3,3),dpra(3,3),eigen(3)
110C----------------------------------------------------------------
111 !=======================================================================
112 ! - RECOVERING MATERIAL PARAMETERS
113 !=======================================================================
114
115 norder = mat_param%IPARAM(1) ! number of Ogden terms
116 iform = mat_param%IPARAM(3) ! Flag for strain energy density formulation
117!
118 ! Shear hyperelastic modulus parameters and material exponents
119 DO i=1,norder
120 mu(i) = mat_param%UPARAM(i)
121 al(i) = mat_param%UPARAM(10+i)
122 END DO
123 ! Shear modulus
124 gmax = zero
125 DO i=1,norder
126 gmax = gmax + mu(i) * al(i)
127 END DO
128 ! Bulk modulus
129 rbulk = mat_param%UPARAM(21)
130 ! Cutoff stress in tension
131 tensioncut = mat_param%UPARAM(23)
132 ! Bulk function scale factor
133 ffac = mat_param%UPARAM(24)
134 ! Tabulated bulk function ID
135 kfp = ifunc(1)
136c
137 ! Fill strain tensor
138 DO i = 1,nel
139 av(1,i) = epsxx(i)
140 av(2,i) = epsyy(i)
141 av(3,i) = epszz(i)
142 av(4,i) = epsxy(i)*half
143 av(5,i) = epsyz(i)*half
144 av(6,i) = epszx(i)*half
145 ENDDO
146c
147 ! Compute principal strains
148 ! Note: in simple precision, principal strains are computed
149 ! with double precision
150 IF (iresp == 1) THEN
151 CALL valpvecdp(av,evv,dirprv,nel)
152 ELSE
153 CALL valpvec(av,evv,dirprv,nel)
154 ENDIF
155c
156 ! Compute principal stretches depending on strain formulation
157 ! (Principal stretches = lambda_i)
158 DO i = 1,nel
159 ! -> Logarithmic strains
160 IF (ismstr == 0 .OR. ismstr == 2 .OR. ismstr == 4) THEN
161 ev(1,i) = exp(evv(1,i))
162 ev(2,i) = exp(evv(2,i))
163 ev(3,i) = exp(evv(3,i))
164 ! -> Green-Lagrange strains
165 ELSEIF (ismstr == 10 .OR. ismstr == 12) THEN
166 ev(1,i) = sqrt(evv(1,i)+ one)
167 ev(2,i) = sqrt(evv(2,i)+ one)
168 ev(3,i) = sqrt(evv(3,i)+ one)
169 ! -> Engineering strains
170 ELSE
171 ev(1,i) = evv(1,i)+ one
172 ev(2,i) = evv(2,i)+ one
173 ev(3,i) = evv(3,i)+ one
174 ENDIF
175 ENDDO
176c
177 ! Relative volume computation
178 ! (RHO/RHO0) = def(F) with F = Grad(Strain)
179 DO i = 1,nel
180 rv(i) = (ev(1,i)*ev(2,i)*ev(3,i))
181 ENDDO
182c
183 ! Compute normalized (deviatoric) stretches and bulk modulus
184 ! (Deviatoric principal stretches = lambda_bar_i)
185 DO i = 1,nel
186 ! Deviatoric stretches
187 evm(1,i) = ev(1,i)*rv(i)**(-third)
188 evm(2,i) = ev(2,i)*rv(i)**(-third)
189 evm(3,i) = ev(3,i)*rv(i)**(-third)
190c
191 ! Tabulated bulk modulus with respect to relative volume
192 IF (kfp /= 0) THEN
193 p(i) = rbulk*ffac*finter(kfp,rv(i),npf,tf,dpdmu)
194 ! Constant bulk modulus
195 ELSE
196 p(i) = rbulk
197 ENDIF
198 ENDDO
199c
200 !=======================================================================
201 ! - STRESS TENSOR COMPUTATION
202 !=======================================================================
203c
204 ! Isochoric strain energy density derivation
205 ! Note: here, the table DWDL(:,J) corresponds to the product EVM(J)*DWDEVM(:,J)
206 DO i=1,nel
207 dwdl(1,i) = zero
208 dwdl(2,i) = zero
209 dwdl(3,i) = zero
210 DO j=1,norder
211 dwdl(1,i) = dwdl(1,i) + mu(j)*(evm(1,i)**al(j) - one)
212 dwdl(2,i) = dwdl(2,i) + mu(j)*(evm(2,i)**al(j) - one)
213 dwdl(3,i) = dwdl(3,i) + mu(j)*(evm(3,i)**al(j) - one)
214 END DO
215 ENDDO
216c
217 ! Volumic strain energy density derivation + trace of isochoric derivative
218 ! (Volumic strain energy density equation depending on formulation flag)
219 DO i=1,nel
220 ! -> Standard strain energy density
221 IF (iform == 1) THEN
222 dwdrv(i) = p(i)*(rv(i) - one)
223 ! -> Modified strain energy density
224 ELSEIF (iform == 2) THEN
225 dwdrv(i) = p(i)*(one - one/rv(i))
226 ENDIF
227 sumdwdl(i) = (dwdl(1,i) + dwdl(2,i) + dwdl(3,i))*third
228 ENDDO
229c
230 ! Cauchy principal stresses
231 DO i = 1,nel
232 DO j = 1,3
233 sigprv(j,i) = (dwdl(j,i)-(sumdwdl(i)-rv(i)*dwdrv(i)))/rv(i)
234 ENDDO
235 ENDDO
236c
237 ! Biot stress tensor for tension cutoff only
238 DO i=1,nel
239 t1(1,i) = rv(i)*sigprv(1,i)/ev(1,i)
240 t1(2,i) = rv(i)*sigprv(2,i)/ev(2,i)
241 t1(3,i) = rv(i)*sigprv(3,i)/ev(3,i)
242 ENDDO
243c
244 ! Tension cutoff stress
245 DO i=1,nel
246 DO j=1,3
247 IF (off(i) == zero .OR. t1(j,i) > abs(tensioncut)) THEN
248 sigprv(1,i) = zero
249 sigprv(2,i) = zero
250 sigprv(3,i) = zero
251 off(i) = zero
252 ENDIF
253 ENDDO
254 ENDDO
255c
256 ! Stress tensor, soundspeed and user-variables
257 DO i=1,nel
258c
259 ! Transform principal Cauchy stresses to Global directions
260 signxx(i) = dirprv(1,1,i)*dirprv(1,1,i)*sigprv(1,i)
261 . + dirprv(1,2,i)*dirprv(1,2,i)*sigprv(2,i)
262 . + dirprv(1,3,i)*dirprv(1,3,i)*sigprv(3,i)
263c
264 signyy(i) = dirprv(2,2,i)*dirprv(2,2,i)*sigprv(2,i)
265 . + dirprv(2,3,i)*dirprv(2,3,i)*sigprv(3,i)
266 . + dirprv(2,1,i)*dirprv(2,1,i)*sigprv(1,i)
267c
268 signzz(i) = dirprv(3,3,i)*dirprv(3,3,i)*sigprv(3,i)
269 . + dirprv(3,1,i)*dirprv(3,1,i)*sigprv(1,i)
270 . + dirprv(3,2,i)*dirprv(3,2,i)*sigprv(2,i)
271c
272 signxy(i) = dirprv(1,1,i)*dirprv(2,1,i)*sigprv(1,i)
273 . + dirprv(1,2,i)*dirprv(2,2,i)*sigprv(2,i)
274 . + dirprv(1,3,i)*dirprv(2,3,i)*sigprv(3,i)
275c
276 signyz(i) = dirprv(2,2,i)*dirprv(3,2,i)*sigprv(2,i)
277 . + dirprv(2,3,i)*dirprv(3,3,i)*sigprv(3,i)
278 . + dirprv(2,1,i)*dirprv(3,1,i)*sigprv(1,i)
279c
280 signzx(i) = dirprv(3,3,i)*dirprv(1,3,i)*sigprv(3,i)
281 . + dirprv(3,1,i)*dirprv(1,1,i)*sigprv(1,i)
282 . + dirprv(3,2,i)*dirprv(1,2,i)*sigprv(2,i)
283c
284 ! Save user variables for post-processing
285 uvar(i,1) = max(sigprv(1,i),sigprv(2,i),sigprv(3,i))
286 uvar(i,2) = min(sigprv(1,i),sigprv(2,i),sigprv(3,i))
287 uvar(i,3) = off(i)
288 ! Soundspeed computation
289 ! Note: for Iform = 2, the constant bulk is assumed for the soundspeed computation
290 ! (if stability issue is observed, switch to non-linear bulk)
291 soundsp(i) = sqrt((two_third*gmax+p(i))/rho(i))
292 ! Viscosity parameter set to zero
293 viscmax(i) = zero
294 ENDDO
295 END
296
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
valpvecdp
subroutine valpvecdp(sig, val, vec, nel)
Definition matrix.F:533
valpvec
subroutine valpvec(sig, val, vec, nel)
Definition matrix.F:215
sigeps42
subroutine sigeps42(mat_param, nel, nuvar, nfunc, ifunc, npf, tf, time, timestep, rho0, rho, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, depsxx, depsyy, depszz, depsxy, depsyz, depszx, epsxx, epsyy, epszz, epsxy, epsyz, epszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, soundsp, viscmax, uvar, off, wxxdt, wyydt, wzzdt, ismstr, mfxx, mfxy, mfxz, mfyx, mfyy, mfyz, mfzx, mfzy, mfzz)
Definition sigeps42.F:44