OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
visc_prony.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!|| visc_prony ../engine/source/materials/visc/visc_prony.F
25!||--- called by ------------------------------------------------------
26!|| viscmain ../engine/source/materials/visc/viscmain.F
27!||--- uses -----------------------------------------------------
28!|| visc_param_mod ../common_source/modules/mat_elem/visc_param_mod.F90
29!||====================================================================
30 SUBROUTINE visc_prony(VISC ,NPRONY ,NEL ,NVARVIS ,UVARVIS ,
31 . EPSPXX ,EPSPYY ,EPSPZZ ,EPSPXY ,EPSPYZ ,EPSPZX ,
32 . SV1 ,SV2 ,SV3 ,SV4 ,SV5 ,SV6 ,
33 . TIMESTEP,RHO ,VISCMAX ,SOUNDSP ,NVAR_DAMP)
34C-----------------------------------------------
35C M o d u l e s
36C-----------------------------------------------
37 USE visc_param_mod
38C-----------------------------------------------
39C I m p l i c i t T y p e s
40C-----------------------------------------------
41#include "implicit_f.inc"
42C---------+---------+---+---+--------------------------------------------
43C I N P U T A r g u m e n t s
44C-----------------------------------------------
45 INTEGER ,INTENT(IN) :: NEL,NVARVIS,NPRONY
46 INTEGER ,INTENT(IN) :: NVAR_DAMP
47 my_real
48 . timestep
49 my_real
50 . rho(nel),
51 . epspxx(nel),epspyy(nel),epspzz(nel),
52 . epspxy(nel),epspyz(nel),epspzx(nel),
53 . sv1(nel),sv2(nel),sv3(nel),sv4(nel),sv5(nel),sv6(nel)
54 TYPE(visc_param_) ,INTENT(IN) :: VISC
55C-----------------------------------------------
56C O U T P U T A r g u m e n t s
57C-----------------------------------------------
58 my_real
59 . viscmax(nel),soundsp(nel)
60C-----------------------------------------------
61C I N P U T O U T P U T A r g u m e n t s
62C-----------------------------------------------
63 my_real ,INTENT(INOUT) :: uvarvis(nel,nvarvis)
64C-----------------------------------------------
65C L o c a l V a r i a b l e s
66C-----------------------------------------------
67 INTEGER I,J,II,IMOD
68 my_real
69 . KV0,G,DAV
70 my_real, DIMENSION(NEL) :: p,epxx,epyy,epzz,trace
71 my_real
72 . aa(nprony),bb(nprony),gv(nprony),beta(nprony),h0(6),h(6),
73 . aak(nprony),bbk(nprony),betak(nprony),kv(nprony),hp0,
74 . rbulk,hp
75C=======================================================================
76 g = zero
77 rbulk = zero
78 imod = visc%IPARAM(2)
79 kv0 = visc%UPARAM(1)
80 viscmax(1:nel) = zero
81
82 DO j=1,nprony
83 gv(j) = visc%UPARAM(1 + j)
84 beta(j) = visc%UPARAM(1 + j + nprony )
85 kv(j) = visc%UPARAM(1 + j + nprony*2)
86 betak(j) = visc%UPARAM(1 + j + nprony*3)
87 g = g + gv(j)
88 rbulk = rbulk + kv(j)
89 aa(j) = exp(-beta(j)*timestep)
90 bb(j) = two*timestep*gv(j)*exp(-half*beta(j)*timestep)
91C
92 aak(j) = exp(-betak(j)*timestep)
93 bbk(j) = timestep*kv(j)*exp(-half*betak(j)*timestep)
94 ENDDO
95C
96 IF(imod > 0) THEN
97 DO i=1,nel
98c spheric part
99 trace(i) = -(epspxx(i) + epspyy(i) + epspzz(i))
100 dav = third*trace(i)
101 p(i) = zero
102
103c deviatoric part
104 epxx(i) = epspxx(i) + dav
105 epyy(i) = epspyy(i) + dav
106 epzz(i) = epspzz(i) + dav
107C
108 sv1(i) = zero
109 sv2(i) = zero
110 sv3(i) = zero
111 sv4(i) = zero
112 sv5(i) = zero
113 sv6(i) = zero
114 ENDDO
115C
116 DO j= 1,nprony
117 ii = (j-1)*(nvarvis-nvar_damp)/nprony
118 DO i=1,nel
119C
120 h0(1) = uvarvis(i,ii + 1)
121 h0(2) = uvarvis(i,ii + 2)
122 h0(3) = uvarvis(i,ii + 3)
123 h0(4) = uvarvis(i,ii + 4)
124 h0(5) = uvarvis(i,ii + 5)
125 h0(6) = uvarvis(i,ii + 6)
126 hp0 = uvarvis(i,ii + 7)
127C
128 h(1) = aa(j)*h0(1) + bb(j)*epxx(i)
129 h(2) = aa(j)*h0(2) + bb(j)*epyy(i)
130 h(3) = aa(j)*h0(3) + bb(j)*epzz(i)
131 h(4) = aa(j)*h0(4) + half*bb(j)*epspxy(i)
132 h(5) = aa(j)*h0(5) + half*bb(j)*epspyz(i)
133 h(6) = aa(j)*h0(6) + half*bb(j)*epspzx(i)
134 hp = aak(j)*hp0 + bbk(j)*trace(i)
135C
136 uvarvis(i,ii + 1) = h(1)
137 uvarvis(i,ii + 2) = h(2)
138 uvarvis(i,ii + 3) = h(3)
139 uvarvis(i,ii + 4) = h(4)
140 uvarvis(i,ii + 5) = h(5)
141 uvarvis(i,ii + 6) = h(6)
142 uvarvis(i,ii + 7) = hp
143c
144 sv1(i) = sv1(i) + h(1)
145 sv2(i) = sv2(i) + h(2)
146 sv3(i) = sv3(i) + h(3)
147 sv4(i) = sv4(i) + h(4)
148 sv5(i) = sv5(i) + h(5)
149 sv6(i) = sv6(i) + h(6)
150 p(i) = p(i) + hp
151 ENDDO
152 ENDDO
153c
154 DO i=1,nel
155 sv1(i) = sv1(i) - p(i)
156 sv2(i) = sv2(i) - p(i)
157 sv3(i) = sv3(i) - p(i)
158C
159 soundsp(i) = sqrt(soundsp(i)**2 + (four_over_3*g + rbulk)/rho(i))
160 ENDDO ! I=1,NEL
161 ELSE
162 DO i=1,nel
163c spheric part
164 trace(i) = epspxx(i) + epspyy(i) + epspzz(i)
165 dav = third*trace(i)
166 p(i) = -kv0*trace(i)
167c deviatoric part
168 epxx(i) = epspxx(i) - dav
169 epyy(i) = epspyy(i) - dav
170 epzz(i) = epspzz(i) - dav
171C
172 sv1(i) = zero
173 sv2(i) = zero
174 sv3(i) = zero
175 sv4(i) = zero
176 sv5(i) = zero
177 sv6(i) = zero
178 ENDDO
179C
180 DO j= 1,nprony
181 ii = (j-1)*(nvarvis-nvar_damp)/nprony
182C
183 DO i=1,nel
184 h0(1) = uvarvis(i,ii + 1)
185 h0(2) = uvarvis(i,ii + 2)
186 h0(3) = uvarvis(i,ii + 3)
187 h0(4) = uvarvis(i,ii + 4)
188 h0(5) = uvarvis(i,ii + 5)
189 h0(6) = uvarvis(i,ii + 6)
190 hp0 = uvarvis(i,ii + 7)
191C
192 h(1) = aa(j)*h0(1) + bb(j)*epxx(i)
193 h(2) = aa(j)*h0(2) + bb(j)*epyy(i)
194 h(3) = aa(j)*h0(3) + bb(j)*epzz(i)
195 h(4) = aa(j)*h0(4) + half*bb(j)*epspxy(i)
196 h(5) = aa(j)*h0(5) + half*bb(j)*epspyz(i)
197 h(6) = aa(j)*h0(6) + half*bb(j)*epspzx(i)
198C
199 uvarvis(i,ii + 1) = h(1)
200 uvarvis(i,ii + 2) = h(2)
201 uvarvis(i,ii + 3) = h(3)
202 uvarvis(i,ii + 4) = h(4)
203 uvarvis(i,ii + 5) = h(5)
204 uvarvis(i,ii + 6) = h(6)
205c
206 sv1(i) = sv1(i) + h(1)
207 sv2(i) = sv2(i) + h(2)
208 sv3(i) = sv3(i) + h(3)
209 sv4(i) = sv4(i) + h(4)
210 sv5(i) = sv5(i) + h(5)
211 sv6(i) = sv6(i) + h(6)
212 ENDDO
213 ENDDO
214c
215 DO i=1,nel
216 sv1(i) = sv1(i) - p(i)
217 sv2(i) = sv2(i) - p(i)
218 sv3(i) = sv3(i) - p(i)
219C
220 soundsp(i) = sqrt(soundsp(i)**2 + (four_over_3*g + rbulk)/rho(i))
221 ENDDO ! I=1,NEL
222 ENDIF
223c------------
224 RETURN
225 END
visc_prony
subroutine visc_prony(visc, nprony, nel, nvarvis, uvarvis, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, sv1, sv2, sv3, sv4, sv5, sv6, timestep, rho, viscmax, soundsp, nvar_damp)
Definition visc_prony.F:34