OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
amomt2.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!|| amomt2 ../engine/source/ale/ale2d/amomt2.F
25!||--- called by ------------------------------------------------------
26!|| qforc2 ../engine/source/elements/solid_2d/quad/qforc2.F
27!||--- calls -----------------------------------------------------
28!|| upwind ../engine/source/elements/solid/solide/upwind.F
29!||--- uses -----------------------------------------------------
30!|| ale_mod ../common_source/modules/ale/ale_mod.F
31!||====================================================================
32 SUBROUTINE amomt2(
33 . PM , W , RHO,
34 . Y1, Y2, Y3, Y4, Z1, Z2, Z3, Z4,
35 . T11, T12 , T13, T14, T21 , T22, T23 , T24,
36 . PY1, PY2 , PZ1, PZ2, AIREM,
37 . VY1, VY2 , VY3, VY4, VZ1 , VZ2, VZ3 , VZ4,
38 . DYY, DZZ , DYZ, DZY,
39 . NC1, NC2 , NC3, NC4, MAT , OFF,
40 . QMV, BUFMAT, DELTAX, VIS, IPM)
41C-----------------------------------------------
42C D e s c r i p t i o n
43C-----------------------------------------------
44C This subroutines computes convective term
45C from momentum continuity equation in 2D.
46C it results into a forces F1,F2 computed at cell
47C centroid.
48C-----------------------------------------------
49C M o d u l e s
50C-----------------------------------------------
51 USE ale_mod
52C-----------------------------------------------
53C I m p l i c i t T y p e s
54C-----------------------------------------------
55#include "implicit_f.inc"
56C-----------------------------------------------
57C G l o b a l P a r a m e t e r s
58C-----------------------------------------------
59#include "mvsiz_p.inc"
60C-----------------------------------------------
61C C o m m o n B l o c k s
62C-----------------------------------------------
63#include "com04_c.inc"
64#include "vect01_c.inc"
65#include "param_c.inc"
66C-----------------------------------------------
67C D u m m y A r g u m e n t s
68C-----------------------------------------------
69 my_real
70 . pm(npropm,nummat), w(3,numnod), rho(*),
71 . y1(*), y2(*), y3(*), y4(*), z1(*), z2(*), z3(*), z4(*),
72 . t11(*), t12(*), t13(*), t14(*), t21(*), t22(*), t23(*),t24(*),
73 . py1(*), py2(*), pz1(*), pz2(*), airem(*),
74 . vy1(*), vy2(*), vy3(*), vy4(*),
75 . vz1(*), vz2(*), vz3(*), vz4(*),
76 . dyy(*), dzz(*), dyz(*), dzy(*), off(*),qmv(8,*),bufmat(*), deltax(*), vis(*)
77 INTEGER MAT(*), NC1(*), NC2(*), NC3(*), NC4(*)
78 INTEGER,INTENT(IN) :: IPM(NPROPMI,NUMMAT)
79C-----------------------------------------------
80C L o c a l V a r i a b l e s
81C-----------------------------------------------
82 INTEGER :: I,IADBUF,IFLG
83 my_real :: F1(MVSIZ), F2(MVSIZ),A1(MVSIZ), A2(MVSIZ), XMS(MVSIZ), GAMMA(MVSIZ),VDY(MVSIZ), VDZ(MVSIZ)
84 my_real :: FAC,AAA,DM,DMY,DMZ
85 my_real :: fvy1(mvsiz), fvz1(mvsiz), fvy2(mvsiz), fvz2(mvsiz), fvy3(mvsiz), fvz3(mvsiz), fvy4(mvsiz), fvz4(mvsiz)
86 my_real :: pc, pd, ps, pp
87 my_real
88 . y12(mvsiz), z12(mvsiz), y23(mvsiz), z23(mvsiz),
89 . y34(mvsiz), z34(mvsiz), y41(mvsiz), z41(mvsiz),
90 . y24(mvsiz), z24(mvsiz), y13(mvsiz), z13(mvsiz),
91 . dyyp(mvsiz), dzzp(mvsiz), dyzp(mvsiz), dzyp(mvsiz),
92 . f1p(mvsiz), f2p(mvsiz), a3(mvsiz),
93 . s(3,mvsiz), t(3,mvsiz), tr(mvsiz)
94
95C-----------------------------------------------
96C S o u r c e C o d e
97C-----------------------------------------------
98 !-----------------------------------------------
99 ! RELATIVE VELOCITY
100 !-----------------------------------------------
101 DO i=lft,llt
102 ! Save fluid velocity
103 fvy1(i) = vy1(i)
104 vy1(i) = vy1(i) - w(2,nc1(i))
105 ! Save fluid velocity
106 fvz1(i) = vz1(i)
107 vz1(i) = vz1(i) - w(3,nc1(i))
108 ! Save fluid velocity
109 fvy2(i) = vy2(i)
110 vy2(i) = vy2(i) - w(2,nc2(i))
111 ! Save fluid velocity
112 fvz2(i) = vz2(i)
113 vz2(i) = vz2(i) - w(3,nc2(i))
114 ! Save fluid velocity
115 fvy3(i) = vy3(i)
116 vy3(i) = vy3(i) - w(2,nc3(i))
117 ! Save fluid velocity
118 fvz3(i) = vz3(i)
119 vz3(i) = vz3(i) - w(3,nc3(i))
120 ! Save fluid velocity
121 fvy4(i) = vy4(i)
122 vy4(i) = vy4(i) - w(2,nc4(i))
123 ! Save fluid velocity
124 fvz4(i) = vz4(i)
125 vz4(i) = vz4(i) - w(3,nc4(i))
126 ENDDO
127
128 !-----------------------------------------------
129 ! CENTROID VALUE (MEAN)
130 !-----------------------------------------------
131 DO i=lft,llt
132 xms(i) = fourth*rho(i)*airem(i)
133 gamma(i)= pm(15,mat(i))
134 ENDDO
135 DO i=lft,llt
136 vdy(i) = fourth*(vy1(i)+vy2(i)+vy3(i)+vy4(i))
137 vdz(i) = fourth*(vz1(i)+vz2(i)+vz3(i)+vz4(i))
138 ENDDO
139
140 !-----------------------------------------------
141 ! LAW51 : convective terme computed like in FVM
142 !-----------------------------------------------
143 ! DMi = 0.5*DT1*QMV(i,I) masse entrant par la face i
144 IF(mtn == 51 .AND. ale%UPWIND%UPWM /= 3) THEN
145 iadbuf = ipm(27,mat(1))
146 iflg = nint(bufmat(31+iadbuf-1))
147 IF(iflg > 1)RETURN
148 IF(ale%UPWIND%UPWM == 0.)THEN
149 DO i=lft,llt
150 gamma(i)= pm(15,mat(i))
151 ENDDO
152 ELSE
153 DO i=lft,llt
154 gamma(i)= ale%UPWIND%CUPWM
155 ENDDO
156 ENDIF
157 DO i=lft,llt
158 aaa = qmv(1,i)+qmv(2,i)+qmv(3,i)+qmv(4,i)
159 aaa = fourth * aaa
160 qmv(1,i) = fourth * (qmv(1,i) - aaa)
161 qmv(2,i) = fourth * (qmv(2,i) - aaa)
162 qmv(3,i) = fourth * (qmv(3,i) - aaa)
163 qmv(4,i) = fourth * (qmv(4,i) - aaa)
164 dmy = zero
165 dmz = zero
166 dm = qmv(4,i)+qmv(1,i)
167 dmy = dmy + vy1(i)*dm
168 dmz = dmz + vz1(i)*dm
169 dm = qmv(1,i)+qmv(2,i)
170 dmy = dmy + vy2(i)*dm
171 dmz = dmz + vz2(i)*dm
172 dm = qmv(2,i)+qmv(3,i)
173 dmy = dmy + vy3(i)*dm
174 dmz = dmz + vz3(i)*dm
175 dm = qmv(3,i)+qmv(4,i)
176 dmy = dmy + vy4(i)*dm
177 dmz = dmz + vz4(i)*dm
178 dmy = -fourth * dmy
179 dmz = -fourth * dmz
180 a1(i) = py1(i)*vdy(i)+pz1(i)*vdz(i)
181 a2(i) = py2(i)*vdy(i)+pz2(i)*vdz(i)
182 a1(i) = sign(gamma(i),a1(i))
183 a2(i) = sign(gamma(i),a2(i))
184 t11(i) = t11(i)+(one+a1(i))*dmy
185 t12(i) = t12(i)+(one+a2(i))*dmy
186 t13(i) = t13(i)+(one-a1(i))*dmy
187 t14(i) = t14(i)+(one-a2(i))*dmy
188 t21(i) = t21(i)+(one+a1(i))*dmz
189 t22(i) = t22(i)+(one+a2(i))*dmz
190 t23(i) = t23(i)+(one-a1(i))*dmz
191 t24(i) = t24(i)+(one-a2(i))*dmz
192 ENDDO
193 ELSEIF (mtn /= 51 .AND. ale%UPWIND%UPWM /= 3) THEN
194 !-----------------------------------------------
195 ! other material laws
196 !-----------------------------------------------
197 DO i=lft,llt
198 f1(i) = (vdy(i)*dyy(i)+vdz(i)*dyz(i))*xms(i)*off(i)
199 f2(i) = (vdy(i)*dzy(i)+vdz(i)*dzz(i))*xms(i)*off(i)
200 ENDDO
201 DO i=lft,llt
202 a1(i) = py1(i)*vdy(i)+pz1(i)*vdz(i)
203 a2(i) = py2(i)*vdy(i)+pz2(i)*vdz(i)
204 ENDDO
205 DO i=lft,llt
206 a1(i) = sign(gamma(i),a1(i))
207 a2(i) = sign(gamma(i),a2(i))
208 ENDDO
209 DO i=lft,llt
210 t11(i) = t11(i)+(one+a1(i))*f1(i)
211 t12(i) = t12(i)+(one+a2(i))*f1(i)
212 t13(i) = t13(i)+(one-a1(i))*f1(i)
213 t14(i) = t14(i)+(one-a2(i))*f1(i)
214
215 t21(i) = t21(i)+(one+a1(i))*f2(i)
216 t22(i) = t22(i)+(one+a2(i))*f2(i)
217 t23(i) = t23(i)+(one-a1(i))*f2(i)
218 t24(i) = t24(i)+(one-a2(i))*f2(i)
219 ENDDO
220
221 ELSEIF (ale%UPWIND%UPWM == 3) THEN
222 pc=fourth
223 pd=one_over_16
224 ps=one_over_16
225 pp=one_over_16
226 DO i=lft,llt
227 y13(i)=y3(i)-y1(i)
228 z13(i)=z3(i)-z1(i)
229 y24(i)=y4(i)-y2(i)
230 z24(i)=z4(i)-z2(i)
231 y12(i)=y2(i)-y1(i)
232 z12(i)=z2(i)-z1(i)
233 y23(i)=y3(i)-y2(i)
234 z23(i)=z3(i)-z2(i)
235 y34(i)=y4(i)-y3(i)
236 z34(i)=z4(i)-z3(i)
237 y41(i)=y1(i)-y4(i)
238 z41(i)=z1(i)-z4(i)
239 ENDDO
240 DO i=lft,llt
241 s(2,i)=half*(y12(i)-y34(i))
242 s(3,i)=half*(z12(i)-z34(i))
243 t(2,i)=half*(-y41(i)+y23(i))
244 t(3,i)=half*(-z41(i)+z23(i))
245 ENDDO
246 !-----------------------------------------------
247 ! SUPG implementation for 2D ALE
248 ! NODE 1
249 !-----------------------------------------------
250 CALL upwind(
251 1 rho, vis, vy1, vy1,
252 2 vz1, s, s, t,
253 3 deltax, gamma, llt)
254 DO i=lft,llt
255 tr(i)=(y41(i)*z12(i)-y12(i)*z41(i))
256 dyyp(i)=(-z24(i)*fvy1(i)-z41(i)*fvy2(i)-z12(i)*fvy4(i))
257 dyzp(i)=( y24(i)*fvy1(i)+y41(i)*fvy2(i)+y12(i)*fvy4(i))
258 dzzp(i)=( y24(i)*fvz1(i)+y41(i)*fvz2(i)+y12(i)*fvz4(i))
259 dzyp(i)=(-z24(i)*fvz1(i)-z41(i)*fvz2(i)-z12(i)*fvz4(i))
260 f1p(i)=rho(i)*(vy1(i)*dyyp(i)+vz1(i)*dyzp(i))
261 f2p(i)=rho(i)*(vy1(i)*dzyp(i)+vz1(i)*dzzp(i))
262 fac=zero
263 IF(tr(i)/=zero)fac=pc*gamma(i)/tr(i)
264 a1(i)=fac*(-z24(i)*vy1(i)+y24(i)*vz1(i))
265 a2(i)=fac*(-z41(i)*vy1(i)+y41(i)*vz1(i))
266 a3(i)=fac*(-z12(i)*vy1(i)+y12(i)*vz1(i))
267
268 t11(i) = t11(i)+(pp+a1(i))*f1p(i)
269 t12(i) = t12(i)+(ps+a2(i))*f1p(i)
270 t13(i) = t13(i)+pd*f1p(i)
271 t14(i) = t14(i)+(ps+a3(i))*f1p(i)
272 t21(i) = t21(i)+(pp+a1(i))*f2p(i)
273 t22(i) = t22(i)+(ps+a2(i))*f2p(i)
274 t23(i) = t23(i)+pd*f2p(i)
275 t24(i) = t24(i)+(ps+a3(i))*f2p(i)
276 ENDDO
277 !-----------------------------------------------
278 ! SUPG implementation for 2D ALE
279 ! NODE 2
280 !-----------------------------------------------
281 CALL upwind(
282 1 rho, vis, vy2, vy2,
283 2 vz2, s, s, t,
284 3 deltax, gamma, llt)
285 DO i=lft,llt
286 tr(i)=(y12(i)*z23(i)-y23(i)*z12(i))
287 dyyp(i)=(-z23(i)*fvy1(i)+z13(i)*fvy2(i)-z12(i)*fvy3(i))
288 dyzp(i)=( y23(i)*fvy1(i)-y13(i)*fvy2(i)+y12(i)*fvy3(i))
289 dzzp(i)=( y23(i)*fvz1(i)-y13(i)*fvz2(i)+y12(i)*fvz3(i))
290 dzyp(i)=(-z23(i)*fvz1(i)+z13(i)*fvz2(i)-z12(i)*fvz3(i))
291 f1p(i)=rho(i)*(vy2(i)*dyyp(i)+vz2(i)*dyzp(i))
292 f2p(i)=rho(i)*(vy2(i)*dzyp(i)+vz2(i)*dzzp(i))
293 fac=zero
294 IF(tr(i)/=zero)fac=pc*gamma(i)/tr(i)
295 a1(i)=fac*(-z23(i)*vy2(i)+y23(i)*vz2(i))
296 a2(i)=fac*( z13(i)*vy2(i)-y13(i)*vz2(i))
297 a3(i)=fac*(-z12(i)*vy2(i)+y12(i)*vz2(i))
298 t11(i) = t11(i)+(ps+a1(i))*f1p(i)
299 t12(i) = t12(i)+(pp+a2(i))*f1p(i)
300 t13(i) = t13(i)+(ps+a3(i))*f1p(i)
301 t14(i) = t14(i)+pd*f1p(i)
302 t21(i) = t21(i)+(ps+a1(i))*f2p(i)
303 t22(i) = t22(i)+(pp+a2(i))*f2p(i)
304 t23(i) = t23(i)+(ps+a3(i))*f2p(i)
305 t24(i) = t24(i)+pd*f2p(i)
306 ENDDO
307 !-----------------------------------------------
308 ! SUPG implementation for 2D ALE
309 ! NODE 3
310 !-----------------------------------------------
311 CALL upwind(
312 1 rho, vis, vy3, vy3,
313 2 vz3, s, s, t,
314 3 deltax, gamma, llt)
315 DO i=lft,llt
316 tr(i)=(y23(i)*z34(i)-y34(i)*z23(i))
317 dyyp(i)=(-z34(i)*fvy2(i)+z24(i)*fvy3(i)-z23(i)*fvy4(i))
318 dyzp(i)=( y34(i)*fvy2(i)-y24(i)*fvy3(i)+y23(i)*fvy4(i))
319 dzzp(i)=( y34(i)*fvz2(i)-y24(i)*fvz3(i)+y23(i)*fvz4(i))
320 dzyp(i)=(-z34(i)*fvz2(i)+z24(i)*fvz3(i)-z23(i)*fvz4(i))
321 f1p(i)=rho(i)*(vy3(i)*dyyp(i)+vz3(i)*dyzp(i))
322 f2p(i)=rho(i)*(vy3(i)*dzyp(i)+vz3(i)*dzzp(i))
323 fac=zero
324 IF(tr(i)/=zero)fac=pc*gamma(i)/tr(i)
325 a1(i)=fac*(-z34(i)*vy3(i)+y34(i)*vz3(i))
326 a2(i)=fac*( z24(i)*vy3(i)-y24(i)*vz3(i))
327 a3(i)=fac*(-z23(i)*vy3(i)+y23(i)*vz3(i))
328 t11(i) = t11(i)+pd*f1p(i)
329 t12(i) = t12(i)+(ps+a1(i))*f1p(i)
330 t13(i) = t13(i)+(pp+a2(i))*f1p(i)
331 t14(i) = t14(i)+(ps+a3(i))*f1p(i)
332 t21(i) = t21(i)+pd*f2p(i)
333 t22(i) = t22(i)+(ps+a1(i))*f2p(i)
334 t23(i) = t23(i)+(pp+a2(i))*f2p(i)
335 t24(i) = t24(i)+(ps+a3(i))*f2p(i)
336 ENDDO
337 !-----------------------------------------------
338 ! SUPG implementation for 2D ALE
339 ! NODE 4
340 !-----------------------------------------------
341 CALL upwind(
342 1 rho, vis, vy4, vy4,
343 2 vz4, s, s, t,
344 3 deltax, gamma, llt)
345 DO i=lft,llt
346 tr(i)=(y34(i)*z41(i)-y41(i)*z34(i))
347 dyyp(i)=(-z34(i)*fvy1(i)-z41(i)*fvy3(i)-z13(i)*fvy4(i))
348 dyzp(i)=( y34(i)*fvy1(i)+y41(i)*fvy3(i)+y13(i)*fvy4(i))
349 dzzp(i)=( y34(i)*fvz1(i)+y41(i)*fvz3(i)+y13(i)*fvz4(i))
350 dzyp(i)=(-z34(i)*fvz1(i)-z41(i)*fvz3(i)-z13(i)*fvz4(i))
351 f1p(i)=rho(i)*(vy4(i)*dyyp(i)+vz4(i)*dyzp(i))
352 f2p(i)=rho(i)*(vy4(i)*dzyp(i)+vz4(i)*dzzp(i))
353 fac=zero
354 IF(tr(i)/=zero)fac=pc*gamma(i)/tr(i)
355 a1(i)=fac*(-z34(i)*vy4(i)+y34(i)*vz4(i))
356 a2(i)=fac*(-z41(i)*vy4(i)+y41(i)*vz4(i))
357 a3(i)=fac*(-z13(i)*vy4(i)+y13(i)*vz4(i))
358 t11(i) = t11(i)+(ps+a1(i))*f1p(i)
359 t12(i) = t12(i)+pd*f1p(i)
360 t13(i) = t13(i)+(ps+a2(i))*f1p(i)
361 t14(i) = t14(i)+(pp+a3(i))*f1p(i)
362 t21(i) = t21(i)+(ps+a1(i))*f2p(i)
363 t22(i) = t22(i)+pd*f2p(i)
364 t23(i) = t23(i)+(ps+a2(i))*f2p(i)
365 t24(i) = t24(i)+(pp+a3(i))*f2p(i)
366 ENDDO
367 ENDIF
368C-----------------------------------------------
369 RETURN
370 END
amomt2
subroutine amomt2(pm, w, rho, y1, y2, y3, y4, z1, z2, z3, z4, t11, t12, t13, t14, t21, t22, t23, t24, py1, py2, pz1, pz2, airem, vy1, vy2, vy3, vy4, vz1, vz2, vz3, vz4, dyy, dzz, dyz, dzy, nc1, nc2, nc3, nc4, mat, off, qmv, bufmat, deltax, vis, ipm)
Definition amomt2.F:41
my_real
#define my_real
Definition cppsort.cpp:32
ale_mod
Definition ale_mod.F:149
ale_mod::ale
type(ale_) ale
Definition ale_mod.F:249
upwind
subroutine upwind(rho, vis, vdx, vdy, vdz, r, s, t, deltax, gam, nel)
Definition upwind.F:35