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