OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
alefvm_stress.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!|| alefvm_stress ../engine/source/ale/alefvm/alefvm_stress.f
25!||--- called by ------------------------------------------------------
26!|| sforc3 ../engine/source/elements/solid/solide/sforc3.F
27!||--- uses -----------------------------------------------------
28!|| alefvm_mod ../common_source/modules/ale/alefvm_mod.F
29!||====================================================================
30 SUBROUTINE alefvm_stress (
31 1 IXS, SIG, QVIS,
32 2 N1X, N2X, N3X, N4X, N5X, N6X,
33 3 N1Y, N2Y, N3Y, N4Y, N5Y, N6Y,
34 4 N1Z, N2Z, N3Z, N4Z, N5Z, N6Z,
35 5 X , IPM, RHO, VOL,
36 . NEL, MOM, SSP, SVIS)
37C-----------------------------------------------
38C D e s c r i p t i o n
39C-----------------------------------------------
40C 'alefvm' is related to a collocated scheme (built from FVM and based on Godunov scheme)
41C which was temporarily introduced for experimental option /INTER/TYPE22 (FSI coupling with cut cell method)
42C This cut cell method is not completed, abandoned, and is not an official option.
43C There is no other use for this scheme which is automatically enabled when /INTER/TYPE22 is defined (INT22>0 => IALEFVM=1).
44C
45C This subroutine is treating an uncut cell.
46C-----------------------------------------------
47C M o d u l e s
48C-----------------------------------------------
49 USE alefvm_mod
50C-----------------------------------------------
51C I m p l i c i t T y p e s
52C-----------------------------------------------
53#include "implicit_f.inc"
54C-----------------------------------------------
55C G l o b a l P a r a m e t e r s
56C-----------------------------------------------
57#include "mvsiz_p.inc"
58C-----------------------------------------------
59C C o m m o n B l o c k s
60C-----------------------------------------------
61#include "vect01_c.inc"
62#include "com01_c.inc"
63#include "param_c.inc"
64C-----------------------------------------------
65C D e s c r i p t i o n
66C-----------------------------------------------
67C This subroutines computes internal forces for
68C finite volume scheme (IALEFVM==1)
69C
70C If option is not detected in input file then
71C subroutine is unplugged
72C-----------------------------------------------
73C D u m m y A r g u m e n t s
74C-----------------------------------------------
75 INTEGER :: NEL
76 INTEGER :: IXS(NIXS,*),IPM(NPROPMI,*)
77 my_real :: SIG(NEL,6),QVIS(*),X(3,*),RHO(*), VOL(*)
78 my_real,INTENT(IN) :: MOM(NEL,3), SSP(*)
79 my_real :: N1X(*), N2X(*), N3X(*), N4X(*), N5X(*), N6X(*),
80 . N1Y(*), N2Y(*), N3Y(*), N4Y(*), N5Y(*), N6Y(*),
81 . n1z(*), n2z(*), n3z(*), n4z(*), n5z(*), n6z(*)
82 my_real, DIMENSION(MVSIZ,6), INTENT(INOUT) :: svis
83C-----------------------------------------------
84C L o c a l V a r i a b l e s
85C-----------------------------------------------
86 INTEGER :: I, II, IMAT, IALEFVM_FLG
87 INTEGER :: NC1(MVSIZ),NC2(MVSIZ),NC3(MVSIZ),NC4(MVSIZ),NC5(MVSIZ),NC6(MVSIZ),NC7(MVSIZ),NC8(MVSIZ)
88
89 my_real :: p(mvsiz) ,
90 . s1(mvsiz), s2(mvsiz), s3(mvsiz),
91 . s4(mvsiz), s5(mvsiz), s6(mvsiz),
92 . mass
93 my_real x1(mvsiz), x2(mvsiz), x3(mvsiz), x4(mvsiz), x5(mvsiz), x6(mvsiz), x7(mvsiz), x8(mvsiz),
94 . y1(mvsiz), y2(mvsiz), y3(mvsiz), y4(mvsiz), y5(mvsiz), y6(mvsiz), y7(mvsiz), y8(mvsiz),
95 . z1(mvsiz), z2(mvsiz), z3(mvsiz), z4(mvsiz), z5(mvsiz), z6(mvsiz), z7(mvsiz), z8(mvsiz),
96 . u_n(6), surf(6), norm(6)
97
98 INTEGER :: idbf,idbl
99 LOGICAL :: debug_outp
100
101C-----------------------------------------------
102C P r e - C o n d i t i o n s
103C-----------------------------------------------
104 IF(alefvm_param%IEnabled==0)RETURN
105 imat = ixs(1,1+nft)
106 ialefvm_flg = ipm(251,imat)
107 IF(ialefvm_flg <= 1)RETURN
108C-----------------------------------------------
109C S o u r c e L i n e s
110C-----------------------------------------------
111
112 ! attention : il faut que mmain soit traite pour chacun des groupes
113 ! afin de connaitre les tenseurs voisins (SIG,SVIS,QVIS) !
114 ! l'assemblage est ensuite realis apr s tous les appel de
115 ! ALEMAIN>SFORC3() pour chaque groupe.
116
117 !-------------------------------------------------------------!
118 ! TOTAL TENSOR !
119 !-------------------------------------------------------------!
120 DO i=1,nel
121 s1(i) = sig(i,1) + svis(i,1) - qvis(i)
122 s2(i) = sig(i,2) + svis(i,2) - qvis(i)
123 s3(i) = sig(i,3) + svis(i,3) - qvis(i)
124 s4(i) = sig(i,4) + svis(i,4)
125 s5(i) = sig(i,5) + svis(i,5)
126 s6(i) = sig(i,6) + svis(i,6)
127 ENDDO
128
129 DO i=1,nel
130 p(i) = -third*(s1(i)+s2(i)+s3(i))
131 ENDDO
132
133
134 !-------------------------------------------------------------!
135 ! NORMAL VECTOR FOR ALE !
136 !-------------------------------------------------------------!
137 DO i=1,nel
138 ii = i + nft
139 !---8 local node numbers NC1 TO NC8 for solid element I ---!
140 nc1(i)=ixs(2,ii)
141 nc2(i)=ixs(3,ii)
142 nc3(i)=ixs(4,ii)
143 nc4(i)=ixs(5,ii)
144 nc5(i)=ixs(6,ii)
145 nc6(i)=ixs(7,ii)
146 nc7(i)=ixs(8,ii)
147 nc8(i)=ixs(9,ii)
148 !
149 !---Coordinates of the 8 nodes
150 x1(i)=x(1,nc1(i))
151 y1(i)=x(2,nc1(i))
152 z1(i)=x(3,nc1(i))
153 !
154 x2(i)=x(1,nc2(i))
155 y2(i)=x(2,nc2(i))
156 z2(i)=x(3,nc2(i))
157 !
158 x3(i)=x(1,nc3(i))
159 y3(i)=x(2,nc3(i))
160 z3(i)=x(3,nc3(i))
161 !
162 x4(i)=x(1,nc4(i))
163 y4(i)=x(2,nc4(i))
164 z4(i)=x(3,nc4(i))
165 !
166 x5(i)=x(1,nc5(i))
167 y5(i)=x(2,nc5(i))
168 z5(i)=x(3,nc5(i))
169 !
170 x6(i)=x(1,nc6(i))
171 y6(i)=x(2,nc6(i))
172 z6(i)=x(3,nc6(i))
173 !
174 x7(i)=x(1,nc7(i))
175 y7(i)=x(2,nc7(i))
176 z7(i)=x(3,nc7(i))
177 !
178 x8(i)=x(1,nc8(i))
179 y8(i)=x(2,nc8(i))
180 z8(i)=x(3,nc8(i))
181 ENDDO
182 DO i=1,nel
183 ! Face-1
184 n1x(i)=(y3(i)-y1(i))*(z2(i)-z4(i)) - (z3(i)-z1(i))*(y2(i)-y4(i))
185 n1y(i)=(z3(i)-z1(i))*(x2(i)-x4(i)) - (x3(i)-x1(i))*(z2(i)-z4(i))
186 n1z(i)=(x3(i)-x1(i))*(y2(i)-y4(i)) - (y3(i)-y1(i))*(x2(i)-x4(i))
187 ! Face-2
188 n2x(i)=(y7(i)-y4(i))*(z3(i)-z8(i)) - (z7(i)-z4(i))*(y3(i)-y8(i))
189 n2y(i)=(z7(i)-z4(i))*(x3(i)-x8(i)) - (x7(i)-x4(i))*(z3(i)-z8(i))
190 n2z(i)=(x7(i)-x4(i))*(y3(i)-y8(i)) - (y7(i)-y4(i))*(x3(i)-x8(i))
191 ! Face-3
192 n3x(i)=(y6(i)-y8(i))*(z7(i)-z5(i)) - (z6(i)-z8(i))*(y7(i)-y5(i))
193 n3y(i)=(z6(i)-z8(i))*(x7(i)-x5(i)) - (x6(i)-x8(i))*(z7(i)-z5(i))
194 n3z(i)=(x6(i)-x8(i))*(y7(i)-y5(i)) - (y6(i)-y8(i))*(x7(i)-x5(i))
195 ! Face-4
196 n4x(i)=(y2(i)-y5(i))*(z6(i)-z1(i)) - (z2(i)-z5(i))*(y6(i)-y1(i))
197 n4y(i)=(z2(i)-z5(i))*(x6(i)-x1(i)) - (x2(i)-x5(i))*(z6(i)-z1(i))
198 n4z(i)=(x2(i)-x5(i))*(y6(i)-y1(i)) - (y2(i)-y5(i))*(x6(i)-x1(i))
199 ! Face-5
200 n5x(i)=(y7(i)-y2(i))*(z6(i)-z3(i)) - (z7(i)-z2(i))*(y6(i)-y3(i))
201 n5y(i)=(z7(i)-z2(i))*(x6(i)-x3(i)) - (x7(i)-x2(i))*(z6(i)-z3(i))
202 n5z(i)=(x7(i)-x2(i))*(y6(i)-y3(i)) - (y7(i)-y2(i))*(x6(i)-x3(i))
203 ! Face-6
204 n6x(i)=(y8(i)-y1(i))*(z4(i)-z5(i)) - (z8(i)-z1(i))*(y4(i)-y5(i))
205 n6y(i)=(z8(i)-z1(i))*(x4(i)-x5(i)) - (x8(i)-x1(i))*(z4(i)-z5(i))
206 n6z(i)=(x8(i)-x1(i))*(y4(i)-y5(i)) - (y8(i)-y1(i))*(x4(i)-x5(i))
207 ENDDO
208
209 !--------------------------------------------------------------!
210 ! PREPARING BUFFER FOR ALEFVM_SFINT3 !
211 ! 0 : cell centroid data !
212 ! 1:3 : INTEGRAL ON EACH FACE from Integral(DIV(SIGMA),Volume)!
213 ! 4 : max(u+x,u-c) !
214 !--------------------------------------------------------------!
215 DO i=1, nel
216 ii = i + nft
217 mass = rho(i)*vol(i)
218 norm(1) = sqrt(n1x(i)*n1x(i) + n1y(i)*n1y(i) + n1z(i)*n1z(i))
219 norm(2) = sqrt(n2x(i)*n2x(i) + n2y(i)*n2y(i) + n2z(i)*n2z(i))
220 norm(3) = sqrt(n3x(i)*n3x(i) + n3y(i)*n3y(i) + n3z(i)*n3z(i))
221 norm(4) = sqrt(n4x(i)*n4x(i) + n4y(i)*n4y(i) + n4z(i)*n4z(i))
222 norm(5) = sqrt(n5x(i)*n5x(i) + n5y(i)*n5y(i) + n5z(i)*n5z(i))
223 norm(6) = sqrt(n6x(i)*n6x(i) + n6y(i)*n6y(i) + n6z(i)*n6z(i))
224 u_n(1) = (mom(i,1)*n1x(i) + mom(i,2)*n1y(i) + mom(i,3)*n1z(i)) / (mass*norm(1))
225 u_n(2) = (mom(i,1)*n2x(i) + mom(i,2)*n2y(i) + mom(i,3)*n2z(i)) / (mass*norm(2))
226 u_n(3) = (mom(i,1)*n3x(i) + mom(i,2)*n3y(i) + mom(i,3)*n3z(i)) / (mass*norm(3))
227 u_n(4) = (mom(i,1)*n4x(i) + mom(i,2)*n4y(i) + mom(i,3)*n4z(i)) / (mass*norm(4))
228 u_n(5) = (mom(i,1)*n5x(i) + mom(i,2)*n5y(i) + mom(i,3)*n5z(i)) / (mass*norm(5))
229 u_n(6) = (mom(i,1)*n6x(i) + mom(i,2)*n6y(i) + mom(i,3)*n6z(i)) / (mass*norm(6))
230 surf(1) = half*norm(1)
231 surf(2) = half*norm(2)
232 surf(3) = half*norm(3)
233 surf(4) = half*norm(4)
234 surf(5) = half*norm(5)
235 surf(6) = half*norm(6)
236 !---storage
237 ! index1=1 centroid data
238 alefvm_buffer%F_FACE(1 ,1 ,ii) = rho(i)
239 alefvm_buffer%F_FACE(1 ,2 ,ii) = ssp(i)
240 alefvm_buffer%F_FACE(1 ,3 ,ii) = rho(i)*ssp(i)
241 alefvm_buffer%F_FACE(1 ,4 ,ii) = p(i)
242 alefvm_buffer%F_FACE(1 ,5 ,ii) = sqrt(mom(i,1)*mom(i,1)+mom(i,2)*mom(i,2)+mom(i,3)*mom(i,3))/rho(i)/vol(i) / ssp(i) ! material velocity (rho.U/rho) divided by sound speed is mach number : M=U/ssps
243 alefvm_buffer%F_FACE(1 ,6 ,ii) = zero !available
244 ! index1=2 Surface
245 alefvm_buffer%F_FACE(2 ,1:6,ii) = surf(1:6)
246 ! index1=3 Normal Velocity
247 alefvm_buffer%F_FACE(3 ,1:6,ii) = u_n(1:6)
248 ENDDO
249
250 !DEBUG-OUTPUT---------------!
251 if(alefvm_param%IOUTP_STRESS /= 0)then
252 debug_outp = .false.
253 if(alefvm_param%IOUTP_STRESS>0)then
254 do i=lft,llt
255 ii = nft + i
256 if(ixs(11,ii)==alefvm_param%IOUTP_STRESS)THEN
257 debug_outp = .true.
258 idbf = i
259 idbl = i
260 EXIT
261 endif
262 enddo
263 elseif(alefvm_param%IOUTP_STRESS==-1)then
264 debug_outp=.true.
265 idbf = lft
266 idbl = llt
267 endif
268 if(debug_outp)then
269!#!include "lockon.inc"
270 print *, " |----alefvm_stress.F-----|"
271 print *, " | THREAD INFORMATION |"
272 print *, " |------------------------|"
273 print *, " NCYCLE =", ncycle
274 do i=idbf,idbl
275 ii = nft + i
276 print *, " brique=", ixs(11,nft+i)
277 write(*,fmt='(A24,1A26)') " ",
278 . "#-stress Tensor (P+VIS+Q)#"
279
280 write (*,fmt='(A,3E26.14,A)') " | ", sig(i,1),sig(i,4),sig(i,6), " |"
281 write (*,fmt='(A,3E26.14,A)') " P =| ", sig(i,4),sig(i,2),sig(i,5), " |"
282 write (*,fmt='(A,3E26.14,A)') " |_", sig(i,6),sig(i,5),sig(i,3), "_|"
283 write (*,fmt='(A,3E26.14,A)') " | ", svis(i,1),svis(i,4),svis(i,6)," |"
284 write (*,fmt='(A,3E26.14,A)') " VIS =| ", svis(i,4),svis(i,2),svis(i,5)," |"
285 write (*,fmt='(A,3E26.14,A)') " |_", svis(i,6),svis(i,5),svis(i,3),"_|"
286 write (*,fmt='(A,3E26.14,A)') " | ", s1(i),s4(i),s6(i)," |"
287 write (*,fmt='(A,3E26.14,A)') " SIGMA = P+VIS+Q =| ", s4(i),s2(i),s5(i)," |"
288 write (*,fmt='(A,3E26.14,A)') " |_", s6(i),s5(i),s3(i),"_|"
289 write (*,fmt='(A,2E26.14)') " p = ",-third*sum(sig(i,1:3)),p(i)
290 write (*,fmt='(A,1E26.14)') " q = ",qvis(i)
291 write (*,fmt='(A,1E26.14)') " M = ",alefvm_buffer%F_FACE(1 ,5 ,ii)
292 write (*,fmt='(A,1E26.14)') " rho = ",alefvm_buffer%F_FACE(1 ,1 ,ii)
293 write (*,fmt='(A,1E26.14)') " u = ",alefvm_buffer%F_FACE(1 ,5 ,ii)*ssp(i)
294 write (*,fmt='(A,1E26.14)') " ssp = ",alefvm_buffer%F_FACE(1 ,2 ,ii)
295 write (*,fmt='(A,1E26.14)') " z = ",alefvm_buffer%F_FACE(1 ,3 ,ii)
296 write (*,fmt='(A,3E26.14)') " rho.U = ",mom(i,1:3)
297 write (*,fmt='(A,3E26.14)') " V = ",vol(i)
298 write(*,fmt='(A34,6A26)') " ",
299 . "#-------- face_1 ---------","#-------- face_2 ---------",
300 . "#-------- face_3 ---------","#-------- face_4 ---------",
301 . "#-------- face_5 ---------","#-------- face_6 --------#"
302 write (*,fmt='(A,8E26.14)') " SURF =", alefvm_buffer%F_FACE(1,1:6,ii)
303 write (*,fmt='(A,8E26.14)') " <U,n> =", alefvm_buffer%F_FACE(2,1:6,ii)
304 print *, " "
305 enddo
306!#!include "lockoff.inc"
307 endif
308 endif
309 !-----------------------------------------!
310
311 RETURN
312 END
alefvm_stress
subroutine alefvm_stress(ixs, sig, qvis, n1x, n2x, n3x, n4x, n5x, n6x, n1y, n2y, n3y, n4y, n5y, n6y, n1z, n2z, n3z, n4z, n5z, n6z, x, ipm, rho, vol, nel, mom, ssp, svis)
Definition alefvm_stress.F:37
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
alefvm_mod
Definition alefvm_mod.F:88
alefvm_mod::alefvm_buffer
type(alefvm_buffer_), target alefvm_buffer
Definition alefvm_mod.F:120
alefvm_mod::alefvm_param
type(alefvm_param_), target alefvm_param
Definition alefvm_mod.F:121