OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sigeps37_single_cell.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!|| sigeps37_single_cell ../engine/source/interfaces/int22/sigeps37_single_cell.F
25!||--- called by ------------------------------------------------------
26!|| sinit22_fvm ../engine/source/interfaces/int22/sinit22_fvm.F
27!||--- uses -----------------------------------------------------
28!|| alefvm_mod ../common_source/modules/ale/alefvm_mod.F
29!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
30!|| element_mod ../common_source/modules/elements/element_mod.F90
31!|| initbuf_mod ../engine/share/resol/initbuf.F
32!||====================================================================
33 SUBROUTINE sigeps37_single_cell (
34 1 ELBUF_TAB, IXS, BUFMAT , IPARG, IPM,
35 2 IDLOC , NG , brickID, VOL
36 . )
37C-----------------------------------------------
38C D e s c r i p t i o n
39C-----------------------------------------------
40C Interface Type22 (/INTER/TYPE22) is an FSI coupling method based on cut cell method.
41C This experimental cut cell method is not completed, abandoned, and is not an official option.
42C
43C This subroutine is treating material law 37.
44C !UVAR(I,1) : massic percentage of liquid * global density (rho1*V1/V : it needs to give liquid mass multiplying by element volume in aleconve.F)
45C !UVAR(I,2) : density of gas
46C !UVAR(I,3) : density of liquid
47C !UVAR(I,4) : volumetric fraction of liquid
48C !UVAR(I,5) : volumetric fraction of gas
49C-----------------------------------------------
50C M o d u l e s
51C-----------------------------------------------
52 USE initbuf_mod
53 USE elbufdef_mod
54 USE alefvm_mod
55 use element_mod , only : nixs
56C-----------------------------------------------
57C I m p l i c i t T y p e s
58C-----------------------------------------------
59#include "implicit_f.inc"
60C-----------------------------------------------
61C I n c l u d e s
62C-----------------------------------------------
63#include "param_c.inc"
64#include "com01_c.inc"
65#include "com04_c.inc"
66C-----------------------------------------------
67C I N P U T A r g u m e n t s
68C-----------------------------------------------
69 my_real,INTENT(INOUT) :: vol
70 INTEGER,INTENT(IN) :: IDLOC,IPARG(NPARG,*), NG, IXS(NIXS,NUMELS), brickID, IPM(NPROPMI,*)
71 TYPE(BUF_MAT_) ,POINTER :: MBUF
72 my_real,INTENT(IN),TARGET :: bufmat(*)
73 TYPE(g_bufel_) ,POINTER :: GBUF
74 TYPE(l_bufel_) ,POINTER :: LBUF
75 TYPE(elbuf_struct_), TARGET, DIMENSION(NGROUP) :: ELBUF_TAB
76C-----------------------------------------------
77C O U T P U T A r g u m e n t s
78C-----------------------------------------------
79
80C-----------------------------------------------
81C I N P U T O U T P U T A r g u m e n t s
82C-----------------------------------------------
83
84C-----------------------------------------------
85C L o c a l V a r i a b l e s
86C-----------------------------------------------
87 INTEGER NITER, ITER,MT, LLT_, IADBUF
88 my_real
89 . SSP,C1,R1,PMIN,RHO10,RHO20,
90 . rho1,rho2, p,gam,p0,
91 . tol, mas, mas1, mas2,
92 . error, p1,p2,f1,f2,df11,df12,df21,df22,det,drho1,drho2,
93 . psh, ssp1, ssp2, rho, uvar(5)
94C-----------------------------------------------
95C S o u r c e C o d e
96C-----------------------------------------------
97
98 !APPEL SIGEPS37 POUR EQUILIBRE PRESSION
99
100 !
101 llt_ = iparg(2,ng)
102 mbuf => elbuf_tab(ng)%BUFLY(1)%MAT(1,1,1)
103 gbuf => elbuf_tab(ng)%GBUF
104 lbuf => elbuf_tab(ng)%BUFLY(1)%LBUF(1,1,1)
105
106 !print *, IDLOC,NG,LLT_
107
108 uvar(1) = mbuf%VAR((0*llt_ + idloc))
109 uvar(2) = mbuf%VAR((1*llt_ + idloc))
110 uvar(3) = mbuf%VAR((2*llt_ + idloc))
111 uvar(4) = mbuf%VAR((3*llt_ + idloc))
112 uvar(5) = mbuf%VAR((4*llt_ + idloc))
113
114
115 mt = ixs(1,brickid)
116 iadbuf = ipm(7,mt)
117 rho10 = bufmat(iadbuf-1+11)
118 rho20 = bufmat(iadbuf-1+12)
119 r1 = bufmat(iadbuf-1+06)
120 p0 = bufmat(iadbuf-1+09)
121 gam = bufmat(iadbuf-1+05)
122 c1 = bufmat(iadbuf-1+04)
123 pmin = bufmat(iadbuf-1+08)
124 psh = bufmat(iadbuf-1+16)
125
126 rho = gbuf%RHO(idloc)
127
128 uvar(1) = mbuf%VAR((1-1)*llt_ +idloc)
129 uvar(2) = mbuf%VAR((2-1)*llt_ +idloc)
130 uvar(3) = mbuf%VAR((3-1)*llt_ +idloc)
131 uvar(4) = mbuf%VAR((4-1)*llt_ +idloc)
132 uvar(5) = mbuf%VAR((5-1)*llt_ +idloc)
133
134 ! if(ixs(11,brickID)==1804589)then
135 ! print * ,"link switch"
136 ! print *, " ini_iter :rho1.V1/V UVAR(1) =", IXS(11,brickID), UVAR(1)
137 ! print *, " ini_iter :rho2 UVAR(2) =", IXS(11,brickID), UVAR(2)
138 ! print *, " ini_iter :rho1 UVAR(3) =", IXS(11,brickID), UVAR(3)
139 ! print *, " ini_iter :alp1 UVAR(4) =", IXS(11,brickID), UVAR(4)
140 ! print *, " ini_iter :alp2 UVAR(5) =", IXS(11,brickID), UVAR(5)
141 ! print *, " ini_iter :v1 =", IXS(11,brickID), UVAR(4)*VOL
142 ! print *, " ini_iter :v2 =", IXS(11,brickID), UVAR(5)*VOL
143 ! print *, " ini_iter :v=v1+v2 =", IXS(11,brickID), VOL
144 ! endif
145 !------------------------------------!
146 ! EQUILIBRIUM ITERATIONS (T>0) !
147 !------------------------------------!
148 tol = em10
149 niter = 20
150
151 mas = rho * vol
152 mas1 = uvar(1) * vol
153 mas2 = mas - mas1
154 rho2 = uvar( 2)
155 rho1 = uvar( 3)
156 IF (mas1 / mas < em10) THEN
157 !!Phase 2 Only present
158 uvar( 1) = zero
159 uvar( 4) = zero
160 uvar( 5) = one
161 rho2 = mas / vol
162 uvar( 2) = rho2
163 p = p0 * (rho2/rho20)**gam
164 ELSEIF (mas2 / mas < em10) THEN
165 !!Phase 1 only present
166 rho1 = mas / vol
167 uvar( 1) = rho1
168 uvar( 3) = rho1
169 uvar( 4) = one
170 uvar( 5) = zero
171 p = r1 * rho1 - c1 + p0
172 ELSE
173 error = ep30
174 iter = 1
175 DO WHILE(iter < niter .AND. error > tol)
176 p1 = r1 * rho1 - c1 + p0
177 p2 = p0 * (rho2/rho20)**gam
178 f1 = mas1 / rho1 + mas2 / rho2 - vol
179 f2 = p1 - p2
180 df11 = - mas1 / (rho1 * rho1)
181 df12 = - mas2 / (rho2 * rho2)
182 df21 = r1
183 df22 = - gam * p0 / (rho20**gam) * rho2**(gam - one)
184 det = df11 * df22 - df12 * df21
185 drho1 = (-df22 * f1 + df12 * f2) / det
186 drho2 = (df21 * f1 - df11 * f2) / det
187 drho1 = min(three * rho1, max(drho1, - half * rho1))
188 drho2 = min(three * rho2, max(drho2, - half * rho2))
189 rho1 = rho1 + drho1
190 rho2 = rho2 + drho2
191 error = abs(drho1 / rho1) + abs(drho2 / rho2)
192 iter = iter + 1
193 ENDDO
194 IF (abs(p1 - p2) > 1.d-5) THEN
195 !PRINT*, "P1", P1, "P2", P2
196 ENDIF
197 IF (error > tol) THEN
198 print*, "*** WARNING LAW37, convergence tol. ", error, tol
199 ENDIF
200 p = r1 * rho1 - c1 + p0
201 ENDIF
202 ssp1 = r1 * rho1
203 ssp2 = gam * p0 * (rho2/rho20)**gam
204 ! P2 = P0 * ((RHO2/RHO20)**GAM-ONE)
205 !---output
206 uvar(2) = rho2
207 uvar(3) = rho1
208 uvar(4) = uvar(1)/rho1
209 IF(uvar(4)<em20)uvar(4)=zero
210 uvar(5) = one-uvar(4)
211 IF (ssp1 > zero) THEN
212 ssp1 = uvar(4) / ssp1
213 ELSE
214 ssp1 = zero
215 ENDIF
216 IF (ssp2 > zero) THEN
217 ssp2 = uvar(5) / ssp2
218 ELSE
219 ssp2 = zero
220 ENDIF
221 ssp = ssp1 + ssp2
222 ssp = sqrt(one / ssp / rho)
223 !print *, "SSP :", SSP, SQRT(R1), SQRT(GAM*P/RHO2)
224 p = max(pmin, p) + psh
225 !SOUNDSP(I) = SSP
226 uvar(1) = max(zero, uvar(1))
227 uvar(2) = max(zero, uvar(2))
228 uvar(3) = max(zero, uvar(3))
229 uvar(4) = max(zero, uvar(4))
230 uvar(5) = max(zero, uvar(5))
231
232
233 !STORE IN BUFFER
234 mbuf%VAR((1-1)*llt_ +idloc) = uvar(1)
235 mbuf%VAR((2-1)*llt_ +idloc) = uvar(2)
236 mbuf%VAR((3-1)*llt_ +idloc) = uvar(3)
237 mbuf%VAR((4-1)*llt_ +idloc) = uvar(4)
238 mbuf%VAR((5-1)*llt_ +idloc) = uvar(5)
239
240 ! if(ixs(11,brickID)==1804589)then
241 !print *, " end_iter :rho1.V1/V UVAR(1) =", IXS(11,brickID), UVAR(1)
242 !print *, " end_iter :rho2 UVAR(2) =", IXS(11,brickID), UVAR(2)
243 !print *, " end_iter :rho1 UVAR(3) =", IXS(11,brickID), UVAR(3)
244 !print *, " end_iter :alp1 UVAR(4) =", IXS(11,brickID), UVAR(4)
245 !print *, " end_iter :alp2 UVAR(5) =", IXS(11,brickID), UVAR(5)
246 !print *, " end_iter :v1 =", IXS(11,brickID), UVAR(4)*VOL
247 !print *, " end_iter :v2 =", IXS(11,brickID), UVAR(5)*VOL
248 !print *, " end_iter :v=v1+v2 =", IXS(11,brickID), VOL
249 !print *, " end_iter :c =", IXS(11,brickID), SSP
250 !print *, " end_iter :P =", IXS(11,brickID), P
251 ! endif
252
253 lbuf%SSP(idloc) = ssp ;
254
255 alefvm_buffer%FCELL(5, brickid ) = ssp
256 alefvm_buffer%FCELL(6, brickid ) = p
257
258C-----------------------------------------------
259 RETURN
260 END
261
262
263
264
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
alefvm_mod
Definition alefvm_mod.F:88
alefvm_mod::alefvm_buffer
type(alefvm_buffer_), target alefvm_buffer
Definition alefvm_mod.F:120
initbuf_mod
Definition initbuf.F:160
sigeps37_single_cell
subroutine sigeps37_single_cell(elbuf_tab, ixs, bufmat, iparg, ipm, idloc, ng, brickid, vol)
Definition sigeps37_single_cell.F:37