OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
visc_prony_lstrain.F File Reference
#include "implicit_f.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine visc_prony_lstrain (visc, nprony, nel, nvarvis, uvarvis, timestep, mfxx, mfxy, mfxz, mfyx, mfyy, mfyz, mfzx, mfzy, mfzz, s1, s2, s3, s4, s5, s6)

Function/Subroutine Documentation

◆ visc_prony_lstrain()

subroutine visc_prony_lstrain ( type(visc_param_), intent(in) visc,
integer, intent(in) nprony,
integer, intent(in) nel,
integer, intent(in) nvarvis,
dimension(nel,nvarvis), intent(inout) uvarvis,
timestep,
intent(in) mfxx,
intent(in) mfxy,
intent(in) mfxz,
intent(in) mfyx,
intent(in) mfyy,
intent(in) mfyz,
intent(in) mfzx,
intent(in) mfzy,
intent(in) mfzz,
intent(inout) s1,
intent(inout) s2,
intent(inout) s3,
intent(inout) s4,
intent(inout) s5,
intent(inout) s6 )

Definition at line 30 of file visc_prony_lstrain.F.

35C-----------------------------------------------
36C M o d u l e s
37C-----------------------------------------------
38 USE visc_param_mod
39C----------------------------------------------
40C I m p l i c i t T y p e s
41C-----------------------------------------------
42#include "implicit_f.inc"
43C---------+---------+---+---+--------------------------------------------
44C VAR | SIZE |TYP| RW| DEFINITION
45C---------+---------+---+---+--------------------------------------------
46C NEL | 1 | I | R | SIZE OF THE ELEMENT GROUP NEL
47C NUPARAM | 1 | I | R | SIZE OF THE USER PARAMETER ARRAY
48C NUVAR | 1 | I | R | NUMBER OF USER ELEMENT VARIABLES
49C---------+---------+---+---+--------------------------------------------
50C TIMESTEP| 1 | F | R | CURRENT TIME STEP
51C UPARAM | NUPARAM | F | R | USER MATERIAL PARAMETER ARRAY
52C RHO | NEL | F | R | INITIAL DENSITY
53C EPSPXX | NEL | F | R | STRAIN RATE XX
54C EPSPYY | NEL | F | R | STRAIN RATE YY
55C ... | | | |
56C ... | | | |
57C VISC | NEL*6 | F | W | VISCOUS STRESS
58C ... | | | |
59C VISCMAX | NEL | F | W | MAXIMUM DAMPING MODULUS(NEEDED FOR TIME STEP)
60C---------+---------+---+---+--------------------------------------------
61C UVAR |NEL*NUVAR| F |R/W| USER ELEMENT VARIABLE ARRAY
62C---------+---------+---+---+--------------------------------------------
63C I N P U T A r g u m e n t s
64C-----------------------------------------------
65 INTEGER ,INTENT(IN) :: NEL,NVARVIS,NPRONY
66 my_real :: timestep
67 my_real, DIMENSION(NEL), INTENT(IN) :: mfxx,mfxy,mfxz,mfyx,mfyy,mfyz ,
68 . mfzx,mfzy,mfzz
69 my_real, DIMENSION(NEL), INTENT(INOUT) :: s1 ,s2 ,s3 ,s4 ,s5 ,s6
70 my_real ,INTENT(INOUT) :: uvarvis(nel,nvarvis)
71 TYPE(VISC_PARAM_) ,INTENT(IN) :: VISC
72C-----------------------------------------------
73C L o c a l V a r i a b l e s
74C-----------------------------------------------
75 INTEGER :: I,II,N,FLAG_VISC
76 my_real :: dtinv,gstard,gamadstar,ttt,mfxx1,mfyy1,mfzz1
77 my_real, DIMENSION(NPRONY) :: epsild,alphad,betad,gamad,
78 . taud
79
80 my_real, DIMENSION(NEL) :: detf,invdetf,jdefg, p
81 my_real :: defgi(nel,9),defg(nel,9),strs0(nel,6),strs1(nel,6),
82 . qint(nel,6,nprony),qintld(nel,6,nprony), temp(nel,9),
83 . qbar(nel,6),temp6(nel,6),shyp(nel,6)
84C=======================================================================
85 flag_visc = visc%IPARAM(3)
86 gamadstar = zero
87 dtinv = timestep /max(em20,timestep**2)
88 DO i=1,nprony
89 gamad(i) = visc%UPARAM(i)
90 taud(i) = visc%UPARAM(nprony + i)
91 ENDDO
92 gstard = zero
93 DO i = 1, nprony
94 ttt = taud(i)*dtinv
95 epsild(i) = exp(-timestep/taud(i))
96 alphad(i) = one - ttt*(one - epsild(i))
97 betad(i) = one - alphad(i) - epsild(i)
98 alphad(i) = gamad(i)*alphad(i)
99 betad(i) = gamad(i)*betad(i)
100 gstard = gstard + alphad(i)
101 ENDDO
102 gstard = one - gstard
103C Computing the inverse of F
104 DO i=1,nel
105 ! J=DETF
106 mfxx1 = mfxx(i) + one
107 mfyy1 = mfyy(i) + one
108 mfzz1 = mfzz(i) + one
109 !!
110 shyp(i,1) = s1(i)
111 shyp(i,2) = s2(i)
112 shyp(i,3) = s3(i)
113 shyp(i,4) = s4(i)
114 shyp(i,5) = s5(i)
115 shyp(i,6) = s6(i)
116
117 defg(i,1) = mfxx1
118 defg(i,2) = mfyx(i)
119 defg(i,3) = mfzx(i)
120 defg(i,4) = mfxy(i)
121 defg(i,5) = mfyy1
122 defg(i,6) = mfzy(i)
123 defg(i,7) = mfxz(i)
124 defg(i,8) = mfyz(i)
125 defg(i,9) = mfzz1
126 !!
127 detf(i) = mfxx1*mfyy1*mfzz1 + mfxy(i)*mfyz(i)*mfzx(i) +
128 . mfyx(i)*mfzy(i)*mfxz(i) - mfyy1*mfxz(i)*mfzx(i) -
129 . mfxy(i)*mfyx(i)*mfzz1 - mfyz(i)*mfzy(i)*mfxx1
130 !! F^-1
131 invdetf(i) = one/detf(i)
132
133!! rhelp = INVDETF(I)
134 defgi(i,1) = (defg(i,5)*defg(i,9)-defg(i,6)*defg(i,8))*invdetf(i)
135 defgi(i,2) =-(defg(i,2)*defg(i,9)-defg(i,3)*defg(i,8))*invdetf(i)
136 defgi(i,3) = (defg(i,2)*defg(i,6)-defg(i,3)*defg(i,5))*invdetf(i)
137 defgi(i,4) =-(defg(i,4)*defg(i,9)-defg(i,6)*defg(i,7))*invdetf(i)
138 defgi(i,5) = (defg(i,1)*defg(i,9)-defg(i,3)*defg(i,7))*invdetf(i)
139 defgi(i,6) =-(defg(i,1)*defg(i,6)-defg(i,3)*defg(i,4))*invdetf(i)
140 defgi(i,7) = (defg(i,4)*defg(i,8)-defg(i,5)*defg(i,7))*invdetf(i)
141 defgi(i,8) =-(defg(i,1)*defg(i,8)-defg(i,2)*defg(i,7))*invdetf(i)
142 defgi(i,9) = (defg(i,1)*defg(i,5)-defg(i,2)*defg(i,4))*invdetf(i)
143 !
144 jdefg(i) = detf(i)
145 ENDDO
146 IF(flag_visc == 2 ) THEN
147 DO i=1,nel
148 p(i) = third*(s1(i) + s2(i) +s3(i))
149 shyp(i,1) = s1(i) - p(i)
150 shyp(i,2) = s2(i) - p(i)
151 shyp(i,3) = s3(i) - p(i)
152 ENDDO
153 ENDIF
154
155 DO i = 1,nel
156 strs1(i,1) = shyp(i,1)
157 strs1(i,2) = shyp(i,2)
158 strs1(i,3) = shyp(i,3)
159 strs1(i,4) = shyp(i,4)
160 strs1(i,5) = shyp(i,5)
161 strs1(i,6) = shyp(i,6)
162 !!
163 strs0(i,1) = uvarvis(i,1)
164 strs0(i,2) = uvarvis(i,2)
165 strs0(i,3) = uvarvis(i,3)
166 strs0(i,4) = uvarvis(i,4)
167 strs0(i,5) = uvarvis(i,5)
168 strs0(i,6) = uvarvis(i,6)
169 ENDDO
170 !! From cauchy to PK2
171 DO n=1,nprony
172 ii = 6*n ! 6 + 6*(N - 1)
173 DO i=1,nel
174 qint(i,1,n) = uvarvis(i,ii + 1)
175 qint(i,2,n) = uvarvis(i,ii + 2)
176 qint(i,3,n) = uvarvis(i,ii + 3)
177 qint(i,4,n) = uvarvis(i,ii + 4)
178 qint(i,5,n) = uvarvis(i,ii + 5)
179 qint(i,6,n) = uvarvis(i,ii + 6)
180 !!
181 qintld(i,1,n) = epsild(n)*qint(i,1,n) + betad(n)*strs0(i,1)
182 qintld(i,2,n) = epsild(n)*qint(i,2,n) + betad(n)*strs0(i,2)
183 qintld(i,3,n) = epsild(n)*qint(i,3,n) + betad(n)*strs0(i,3)
184 qintld(i,4,n) = epsild(n)*qint(i,4,n) + betad(n)*strs0(i,4)
185 qintld(i,5,n) = epsild(n)*qint(i,5,n) + betad(n)*strs0(i,5)
186 qintld(i,6,n) = epsild(n)*qint(i,6,n) + betad(n)*strs0(i,6)
187 ENDDO
188 ENDDO
189 DO i=1,nel
190 ! f-1*sig
191 temp(i,1) = defgi(i,1)*strs1(i,1) + defgi(i,4)*strs1(i,4) +
192 . defgi(i,7)*strs1(i,6)
193 temp(i,2) = defgi(i,2)*strs1(i,1) + defgi(i,5)*strs1(i,4) +
194 . defgi(i,8)*strs1(i,6)
195 temp(i,3) = defgi(i,3)*strs1(i,1) + defgi(i,6)*strs1(i,4) +
196 . defgi(i,9)*strs1(i,6)
197 temp(i,4) = defgi(i,1)*strs1(i,4) + defgi(i,4)*strs1(i,2) +
198 . defgi(i,7)*strs1(i,5)
199 temp(i,5) = defgi(i,2)*strs1(i,4) + defgi(i,5)*strs1(i,2) +
200 . defgi(i,8)*strs1(i,5)
201 temp(i,6) = defgi(i,3)*strs1(i,4) + defgi(i,6)*strs1(i,2) +
202 . defgi(i,9)*strs1(i,5)
203 temp(i,7) = defgi(i,1)*strs1(i,6) + defgi(i,4)*strs1(i,5) +
204 . defgi(i,7)*strs1(i,3)
205 temp(i,8) = defgi(i,2)*strs1(i,6) + defgi(i,5)*strs1(i,5) +
206 . defgi(i,8)*strs1(i,3)
207 temp(i,9) = defgi(i,3)*strs1(i,6) + defgi(i,6)*strs1(i,5) +
208 . defgi(i,9)*strs1(i,3)
209c
210 strs0(i,1) = jdefg(i)*(temp(i,1)*defgi(i,1) + temp(i,4)*defgi(i,4)
211 . + temp(i,7)*defgi(i,7))
212 strs0(i,2) = jdefg(i)*(temp(i,2)*defgi(i,2) + temp(i,5)*defgi(i,5)
213 . + temp(i,8)*defgi(i,8))
214 strs0(i,3) = jdefg(i)*(temp(i,3)*defgi(i,3) + temp(i,6)*defgi(i,6)
215 . + temp(i,9)*defgi(i,9))
216 strs0(i,4) = jdefg(i)*(temp(i,2)*defgi(i,1) + temp(i,5)*defgi(i,4)
217 . + temp(i,8)*defgi(i,7))
218 strs0(i,5) = jdefg(i)*(temp(i,3)*defgi(i,2) + temp(i,6)*defgi(i,5)
219 . + temp(i,9)*defgi(i,8))
220 strs0(i,6) = jdefg(i)*(temp(i,3)*defgi(i,1) + temp(i,6)*defgi(i,4)
221 . + temp(i,9)*defgi(i,7))
222 uvarvis(i,1) = strs0(i,1)
223 uvarvis(i,2) = strs0(i,2)
224 uvarvis(i,3) = strs0(i,3)
225 uvarvis(i,4) = strs0(i,4)
226 uvarvis(i,5) = strs0(i,5)
227 uvarvis(i,6) = strs0(i,6)
228 ENDDO
229 DO n = 1, nprony
230 ii = 6*n ! 6 + 6*(N-1)
231 DO i = 1,nel
232 qint(i,1,n) = qintld(i,1,n) + alphad(n)*strs0(i,1)
233 qint(i,2,n) = qintld(i,2,n) + alphad(n)*strs0(i,2)
234 qint(i,3,n) = qintld(i,3,n) + alphad(n)*strs0(i,3)
235 qint(i,4,n) = qintld(i,4,n) + alphad(n)*strs0(i,4)
236 qint(i,5,n) = qintld(i,5,n) + alphad(n)*strs0(i,5)
237 qint(i,6,n) = qintld(i,6,n) + alphad(n)*strs0(i,6)
238 !!
239 uvarvis(i,ii + 1) = qint(i,1,n)
240 uvarvis(i,ii + 2) = qint(i,2,n)
241 uvarvis(i,ii + 3) = qint(i,3,n)
242 uvarvis(i,ii + 4) = qint(i,4,n)
243 uvarvis(i,ii + 5) = qint(i,5,n)
244 uvarvis(i,ii + 6) = qint(i,6,n)
245 ENDDO
246 ENDDO
247 ! total viscous stress
248 DO i = 1,nel
249 qbar(i,1) = zero
250 qbar(i,2) = zero
251 qbar(i,3) = zero
252 qbar(i,4) = zero
253 qbar(i,5) = zero
254 qbar(i,6) = zero
255 ENDDO
256c
257 DO n = 1,nprony
258 DO i = 1,nel
259 temp(i,1) = defg(i,1)*qintld(i,1,n) +
260 . defg(i,4)*qintld(i,4,n) +
261 . defg(i,7)*qintld(i,6,n)
262 temp(i,2) = defg(i,2)*qintld(i,1,n) +
263 . defg(i,5)*qintld(i,4,n) +
264 . defg(i,8)*qintld(i,6,n)
265 temp(i,3) = defg(i,3)*qintld(i,1,n) +
266 . defg(i,6)*qintld(i,4,n) +
267 . defg(i,9)*qintld(i,6,n)
268 temp(i,4) = defg(i,1)*qintld(i,4,n) +
269 . defg(i,4)*qintld(i,2,n) +
270 . defg(i,7)*qintld(i,5,n)
271 temp(i,5) = defg(i,2)*qintld(i,4,n) +
272 . defg(i,5)*qintld(i,2,n) +
273 . defg(i,8)*qintld(i,5,n)
274 temp(i,6) = defg(i,3)*qintld(i,4,n) +
275 . defg(i,6)*qintld(i,2,n) +
276 . defg(i,9)*qintld(i,5,n)
277 temp(i,7) = defg(i,1)*qintld(i,6,n) +
278 . defg(i,4)*qintld(i,5,n) +
279 . defg(i,7)*qintld(i,3,n)
280 temp(i,8) = defg(i,2)*qintld(i,6,n) +
281 . defg(i,5)*qintld(i,5,n) +
282 . defg(i,8)*qintld(i,3,n)
283 temp(i,9) = defg(i,3)*qintld(i,6,n) +
284 . defg(i,6)*qintld(i,5,n) +
285 . defg(i,9)*qintld(i,3,n)
286c
287 temp6(i,1) = temp(i,1)*defg(i,1) +
288 . temp(i,4)*defg(i,4) +
289 . temp(i,7)*defg(i,7)
290 temp6(i,2) = temp(i,2)*defg(i,2) +
291 . temp(i,5)*defg(i,5) +
292 . temp(i,8)*defg(i,8)
293 temp6(i,3) = temp(i,3)*defg(i,3) +
294 . temp(i,6)*defg(i,6) +
295 . temp(i,9)*defg(i,9)
296 temp6(i,4) = temp(i,2)*defg(i,1) +
297 . temp(i,5)*defg(i,4) +
298 . temp(i,8)*defg(i,7)
299 temp6(i,5) = temp(i,3)*defg(i,2) +
300 . temp(i,6)*defg(i,5) +
301 . temp(i,9)*defg(i,8)
302 temp6(i,6) = temp(i,3)*defg(i,1) +
303 . temp(i,6)*defg(i,4) +
304 . temp(i,9)*defg(i,7)
305c
306 qbar(i,1) = qbar(i,1) + temp6(i,1)
307 qbar(i,2) = qbar(i,2) + temp6(i,2)
308 qbar(i,3) = qbar(i,3) + temp6(i,3)
309 qbar(i,4) = qbar(i,4) + temp6(i,4)
310 qbar(i,5) = qbar(i,5) + temp6(i,5)
311 qbar(i,6) = qbar(i,6) + temp6(i,6)
312c
313 ENDDO
314 ENDDO
315c
316c updated cauchy stress
317c
318 DO i = 1,nel
319 strs1(i,1) = gstard*strs1(i,1) - qbar(i,1)/jdefg(i)
320 strs1(i,2) = gstard*strs1(i,2) - qbar(i,2)/jdefg(i)
321 strs1(i,3) = gstard*strs1(i,3) - qbar(i,3)/jdefg(i)
322 strs1(i,4) = gstard*strs1(i,4) - qbar(i,4)/jdefg(i)
323 strs1(i,5) = gstard*strs1(i,5) - qbar(i,5)/jdefg(i)
324 strs1(i,6) = gstard*strs1(i,6) - qbar(i,6)/jdefg(i)
325 !!
326 s1(i) = strs1(i,1)
327 s2(i) = strs1(i,2)
328 s3(i) = strs1(i,3)
329 s4(i) = strs1(i,4)
330 s5(i) = strs1(i,5)
331 s6(i) = strs1(i,6)
332 ENDDO
333 IF(flag_visc == 2) THEN
334 DO i=1,nel
335 s1(i) = s1(i) + p(i)
336 s2(i) = s2(i) + p(i)
337 s3(i) = s3(i) + p(i)
338 ENDDO
339 ENDIF
340c------------
341 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
max
#define max(a, b)
Definition macros.h:21