OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps42.F File Reference
#include "implicit_f.inc"
#include "mvsiz_p.inc"
#include "scr05_c.inc"
#include "tabsiz_c.inc"

Go to the source code of this file.

Functions/Subroutines

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)

Function/Subroutine Documentation

◆ sigeps42()

subroutine sigeps42 ( type(matparam_struct_), intent(inout) mat_param,
integer, intent(in) nel,
integer, intent(in) nuvar,
integer nfunc,
integer, dimension(nfunc) ifunc,
integer, dimension(snpc) npf,
tf,
intent(in) time,
intent(in) timestep,
dimension (nel), intent(in) rho0,
dimension (nel), intent(in) rho,
dimension(nel), intent(in) epspxx,
dimension(nel), intent(in) epspyy,
dimension(nel), intent(in) epspzz,
dimension(nel), intent(in) epspxy,
dimension(nel), intent(in) epspyz,
dimension(nel), intent(in) epspzx,
dimension(nel), intent(in) depsxx,
dimension(nel), intent(in) depsyy,
dimension(nel), intent(in) depszz,
dimension(nel), intent(in) depsxy,
dimension(nel), intent(in) depsyz,
dimension(nel), intent(in) depszx,
dimension (nel), intent(in) epsxx,
dimension (nel), intent(in) epsyy,
dimension (nel), intent(in) epszz,
dimension (nel), intent(in) epsxy,
dimension (nel), intent(in) epsyz,
dimension (nel), intent(in) epszx,
dimension(nel), intent(in) sigoxx,
dimension(nel), intent(in) sigoyy,
dimension(nel), intent(in) sigozz,
dimension(nel), intent(in) sigoxy,
dimension(nel), intent(in) sigoyz,
dimension(nel), intent(in) sigozx,
dimension (nel), intent(out) signxx,
dimension (nel), intent(out) signyy,
dimension(nel), intent(out) signzz,
dimension (nel), intent(out) signxy,
dimension (nel), intent(out) signyz,
dimension(nel), intent(out) signzx,
dimension(nel), intent(out) soundsp,
dimension(nel), intent(out) viscmax,
dimension(nel,nuvar), intent(inout) uvar,
dimension(nel), intent(inout) off,
dimension(mvsiz), intent(in) wxxdt,
dimension(mvsiz), intent(in) wyydt,
dimension(mvsiz), intent(in) wzzdt,
integer, intent(in) ismstr,
dimension(nel), intent(in) mfxx,
dimension(nel), intent(in) mfxy,
dimension(nel), intent(in) mfxz,
dimension(nel), intent(in) mfyx,
dimension(nel), intent(in) mfyy,
dimension(nel), intent(in) mfyz,
dimension(nel), intent(in) mfzx,
dimension(nel), intent(in) mfzy,
dimension(nel), intent(in) mfzz )

Definition at line 33 of file sigeps42.F.

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
my_real
#define my_real
Definition cppsort.cpp:32
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