1!------------------------------------------------------------------------------- 2 3! This file is part of Code_Saturne, a general-purpose CFD tool. 4! 5! Copyright (C) 1998-2021 EDF S.A. 6! 7! This program is free software; you can redistribute it and/or modify it under 8! the terms of the GNU General Public License as published by the Free Software 9! Foundation; either version 2 of the License, or (at your option) any later 10! version. 11! 12! This program is distributed in the hope that it will be useful, but WITHOUT 13! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 14! FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 15! details. 16! 17! You should have received a copy of the GNU General Public License along with 18! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin 19! Street, Fifth Floor, Boston, MA 02110-1301, USA. 20 21!------------------------------------------------------------------------------- 22 23subroutine cfrusb & 24!================ 25 26 ( ifac , & 27 bc_en , bc_pr , bc_vel ) 28 29!=============================================================================== 30! FUNCTION : 31! --------- 32 33! Rusanov flux at the boundary for Euler + Energy 34 35! d rho /dt + div rho u = 0 36! d rho u /dt + div rho u u + grad P = 0 37! d E /dt + div rho u E + div u P = 0 38 39!------------------------------------------------------------------------------- 40! Arguments 41!__________________.____._____.________________________________________________. 42! name !type!mode ! role ! 43!__________________!____!_____!________________________________________________! 44! ifac ! i ! <-- ! face number ! 45!__________________!____!_____!________________________________________________! 46 47! TYPE : E (ENTIER), R (REEL), A (ALPHANUMERIQUE), T (TABLEAU) 48! L (LOGIQUE) .. ET TYPES COMPOSES (EX : TR TABLEAU REEL) 49! MODE : <-- donnee, --> resultat, <-> Donnee modifiee 50! --- tableau de travail 51!=============================================================================== 52 53!=============================================================================== 54! Module files 55!=============================================================================== 56 57use paramx 58use pointe, only:rvoid1 59use numvar 60use optcal 61use cstphy 62use cstnum 63use parall 64use entsor 65use ppppar 66use ppthch 67use ppincl 68use cfpoin 69use mesh 70use field 71use cs_cf_bindings 72 73!=============================================================================== 74 75implicit none 76 77! Arguments 78 79integer ifac 80 81double precision bc_en(nfabor), bc_pr(nfabor), bc_vel(3,nfabor) 82 83! Local variables 84 85integer iel, l_size 86integer ien 87 88double precision inv_surfbn, rnx, rny, rnz 89double precision b_vel_n, c_vel_n, b_masfl, c_masfl, b_c, c_c, b_c2(1), c_c2(1) 90double precision rrus, r_b_masfl, r_b_vel_n 91double precision, dimension(:,:), pointer :: cofacv 92double precision, dimension(:), pointer :: coface 93double precision, dimension(:), pointer :: crom, brom 94double precision, dimension(:,:), pointer :: vel 95double precision, dimension(:), pointer :: cvar_pr, cvar_en, cpro_cp, cpro_cv 96double precision, dimension(:), pointer :: rvoid 97 98double precision c_cp(1), c_cv(1) 99 100!=============================================================================== 101 102rvoid => null() 103 104!=============================================================================== 105! 0. INITIALISATION 106!=============================================================================== 107 108! Map field arrays 109call field_get_val_v(ivarfl(iu), vel) 110 111ien = isca(ienerg) 112 113call field_get_val_s(icrom, crom) 114call field_get_val_s(ibrom, brom) 115 116call field_get_val_s(ivarfl(ipr), cvar_pr) 117call field_get_val_s(ivarfl(ien), cvar_en) 118 119call field_get_coefac_v(ivarfl(iu), cofacv) 120call field_get_coefac_s(ivarfl(ien), coface) 121 122iel = ifabor(ifac) 123 124! Initialize local specific heat values and handle uniform cases 125 126if (icp.ge.0) then 127 call field_get_val_s(icp, cpro_cp) 128 c_cp(1) = cpro_cp(iel) 129else 130 cpro_cp => rvoid1 131 c_cp(1) = 0.d0 132endif 133 134if (icv.ge.0) then 135 call field_get_val_s(icv, cpro_cv) 136 c_cv(1) = cpro_cv(iel) 137else 138 cpro_cv => rvoid1 139 c_cv(1) = 0.d0 140endif 141 142!=============================================================================== 143! 1. COMPUTE VALUES NEEDED FOR RUSANOV SCHEME 144!=============================================================================== 145 146inv_surfbn = 1. / surfbn(ifac) 147rnx = surfbo(1,ifac) * inv_surfbn 148rny = surfbo(2,ifac) * inv_surfbn 149rnz = surfbo(3,ifac) * inv_surfbn 150 151b_vel_n = bc_vel(1,ifac)*rnx + bc_vel(2,ifac)*rny + bc_vel(3,ifac)*rnz 152c_vel_n = vel(1,iel)*rnx + vel(2,iel)*rny + vel(3,iel)*rnz 153b_masfl = brom(ifac)*b_vel_n 154c_masfl = crom(iel)*c_vel_n 155 156l_size = 1 157 158call cs_cf_thermo_c_square(c_cp, c_cv, & 159 bc_pr(ifac:ifac), brom(ifac:ifac), & 160 rvoid, rvoid, rvoid, b_c2, l_size) 161call cs_cf_thermo_c_square(c_cp, c_cv, & 162 cvar_pr(iel:iel), crom(iel:iel), & 163 rvoid, rvoid, rvoid, c_c2, l_size) 164 165b_c = sqrt(b_c2(1)) 166c_c = sqrt(c_c2(1)) 167rrus = max(abs(b_vel_n)+b_c, abs(c_vel_n)+c_c) 168 169! boundary mass flux computed with Rusanov scheme 170r_b_masfl = 0.5d0*(b_masfl+c_masfl) - 0.5d0*rrus*(brom(ifac)-crom(iel)) 171 172! Rusanov normal velocity (computed using boundary density) 173r_b_vel_n = r_b_masfl / brom(ifac) 174 175! Update velocity boundary condition 176bc_vel(1,ifac) = bc_vel(1,ifac) + (r_b_vel_n - b_vel_n)*rnx 177bc_vel(2,ifac) = bc_vel(2,ifac) + (r_b_vel_n - b_vel_n)*rny 178bc_vel(3,ifac) = bc_vel(3,ifac) + (r_b_vel_n - b_vel_n)*rnz 179 180!=============================================================================== 181! 2. CONVECTIVE RUSANOV FLUX 182!=============================================================================== 183 184! Tag the faces where a Rusanov flux is computed 185! The tag will be used in bilsc2 to retrieve the faces where a Rusanov flux 186! has to be imposed 187icvfli(ifac) = 1 188 189! Momentum flux (the centered pressure contribution is directly taken into account 190! in the pressure BC) 191cofacv(1,ifac) = suffbn(ifac)* & 192 0.5d0*( b_masfl*bc_vel(1,ifac) + c_masfl*vel(1,iel) & 193 -rrus*(brom(ifac)*bc_vel(1,ifac) & 194 -crom(iel)*vel(1,iel)) ) 195 196cofacv(2,ifac) = suffbn(ifac)* & 197 0.5d0*( b_masfl*bc_vel(2,ifac) + c_masfl*vel(2,iel) & 198 -rrus*( brom(ifac)*bc_vel(2,ifac) & 199 -crom(iel)*vel(2,iel)) ) 200 201cofacv(3,ifac) = suffbn(ifac)* & 202 0.5d0*( b_masfl*bc_vel(3,ifac) + c_masfl*vel(3,iel) & 203 -rrus*(brom(ifac)*bc_vel(3,ifac) & 204 -crom(iel)*vel(3,iel)) ) 205 206! BC for the pressure gradient in the momentum balance 207bc_pr(ifac) = 0.5d0 * (bc_pr(ifac) + cvar_pr(iel)) 208 209! Total energy flux 210coface(ifac) = suffbn(ifac)* & 211 0.5d0*( b_masfl*bc_en(ifac) + c_masfl*cvar_en(iel) & 212 +b_vel_n*bc_pr(ifac) + c_vel_n*cvar_pr(iel) & 213 -rrus*(brom(ifac)*bc_en(ifac) & 214 -crom(iel)*cvar_en(iel)) ) 215 216return 217 218end subroutine 219