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 23!> \file cavitation.f90 24!> Module for cavitation modeling 25!> 26!> Please refer to the 27!> <a href="../../theory.pdf#cavitation"><b>cavitation model</b></a> section 28!> of the theory guide for more informations. 29 30module cavitation 31 32 !============================================================================= 33 34 use, intrinsic :: iso_c_binding 35 36 implicit none 37 38 !============================================================================= 39 40 !> \addtogroup cavitation 41 !> \{ 42 43 !---------------------------------------------------------------------------- 44 ! Vaporization/condensation model 45 !---------------------------------------------------------------------------- 46 47 !> \addtogroup cav_source_term 48 !> \{ 49 50 !> reference saturation pressure (kg/(m s2)) 51 real(c_double), pointer, save :: presat 52 53 !> reference velocity of the flow (m/s) 54 real(c_double), pointer, save :: uinf 55 56 !> reference length scale of the flow (m) 57 real(c_double), pointer, save :: linf 58 59 !> constant Cdest of the condensation source term (Merkle model) 60 real(c_double), pointer, save :: cdest 61 62 !> constant Cprod of the vaporization source term (Merkle model) 63 real(c_double), pointer, save :: cprod 64 65 !> \} 66 67 !---------------------------------------------------------------------------- 68 ! Interaction with turbulence 69 !---------------------------------------------------------------------------- 70 71 !> \addtogroup cav_turbulence 72 !> \{ 73 74 !> activation of the eddy-viscosity correction (Reboud correction) 75 !> - 1: activated 76 !> - 0: desactivated 77 integer(c_int), pointer, save :: icvevm 78 79 !> constant mcav of the eddy-viscosity correction (Reboud correction) 80 real(c_double), pointer, save :: mcav 81 82 !> \} 83 84 !---------------------------------------------------------------------------- 85 ! Numerical parameters 86 !---------------------------------------------------------------------------- 87 88 !> \addtogroup cav_numerics 89 !> \{ 90 91 !> implicitation in pressure of the vaporization/condensation model 92 !> - 1: activated 93 !> - 0: desactivated 94 integer(c_int), pointer, save :: itscvi 95 96 !> \} 97 98 !> \} 99 100 !============================================================================= 101 102 interface 103 104 !--------------------------------------------------------------------------- 105 106 !> \cond DOXYGEN_SHOULD_SKIP_THIS 107 108 !--------------------------------------------------------------------------- 109 110 ! Interface to C function retrieving pointers to cavitation model indicator 111 ! and parameters 112 113 subroutine cs_f_cavitation_get_pointers(presat, uinf, linf, cdest, cprod, & 114 icvevm, mcav, itscvi) & 115 bind(C, name='cs_f_cavitation_get_pointers') 116 use, intrinsic :: iso_c_binding 117 implicit none 118 type(c_ptr), intent(out) :: presat, uinf, linf, cdest, cprod 119 type(c_ptr), intent(out) :: icvevm, mcav, itscvi 120 end subroutine cs_f_cavitation_get_pointers 121 122 !--------------------------------------------------------------------------- 123 124 !> (DOXYGEN_SHOULD_SKIP_THIS) \endcond 125 126 !--------------------------------------------------------------------------- 127 128 end interface 129 130 !============================================================================= 131 132contains 133 134 !============================================================================= 135 136 !> \brief Initialize Fortran cavitation model API. 137 !> This maps Fortran pointers to global C structure members and indicator. 138 139 subroutine cavitation_model_init 140 141 use, intrinsic :: iso_c_binding 142 143 implicit none 144 145 ! Local variables 146 147 type(c_ptr) :: c_presat, c_uinf, c_linf, c_cdest, c_cprod 148 type(c_ptr) :: c_icvevm, c_mcav, c_itscvi 149 150 call cs_f_cavitation_get_pointers(c_presat, c_uinf, c_linf, & 151 c_cdest, c_cprod, c_icvevm, c_mcav, & 152 c_itscvi) 153 154 call c_f_pointer(c_presat, presat) 155 call c_f_pointer(c_uinf, uinf) 156 call c_f_pointer(c_linf, linf) 157 call c_f_pointer(c_cdest, cdest) 158 call c_f_pointer(c_cprod, cprod) 159 call c_f_pointer(c_icvevm, icvevm) 160 call c_f_pointer(c_mcav, mcav) 161 call c_f_pointer(c_itscvi, itscvi) 162 163 end subroutine cavitation_model_init 164 165 !============================================================================= 166 167 !> \brief Compute the vaporization source term 168 !> \f$ \Gamma_V \left(\alpha, p\right) = m^+ + m^- \f$ using the 169 !> Merkle model: 170 !> \f[ 171 !> m^+ = -\dfrac{C_{prod} \rho_l \min \left( p-p_V,0 \right)\alpha(1-\alpha)} 172 !> {0.5\rho_lu_\infty^2t_\infty}, 173 !> \f] 174 !> \f[ 175 !> m^- = -\dfrac{C_{dest} \rho_v \max \left( p-p_V,0 \right)\alpha(1-\alpha)} 176 !> {0.5\rho_lu_\infty^2t_\infty}, 177 !> \f] 178 !> with \f$ C_{prod}, C_{dest} \f$ empirical constants, 179 !> \f$ t_\infty=l_\infty/u_\infty \f$ a reference time scale and \f$p_V\f$ 180 !> the reference saturation pressure. 181 !> \f$ l_\infty \f$, \f$ u_\infty \f$ and \f$p_V\f$ may be provided by 182 !> the user (user function). 183 !> Note that the r.h.s. of the void fraction transport equation is 184 !> \f$ \Gamma_V/\rho_v \f$. 185 186 !> \param[in] pressure Pressure array 187 !> \param[in] voidf Void fraction array 188 189 subroutine cavitation_compute_source_term(pressure, voidf) 190 191 use optcal 192 use pointe, only: gamcav, dgdpca 193 use mesh, only: ncel, ncelet 194 use vof 195 196 ! Arguments 197 198 double precision pressure(ncelet), voidf(ncelet) 199 200 ! Local variables 201 202 integer iel 203 double precision tinf, cond, cvap, condens, vaporis 204 205 if (iand(ivofmt,VOF_MERKLE_MASS_TRANSFER).ne.0) then 206 207 ! Merkle model 208 209 tinf = linf/uinf 210 211 cond = (cdest*rho2)/(0.5d0*rho1*uinf*uinf*tinf) 212 cvap = (cprod*rho1)/(0.5d0*rho1*uinf*uinf*tinf) 213 214 do iel = 1, ncel 215 condens = -cond*max(0.d0, pressure(iel) - presat) & 216 *voidf(iel)*(1.d0 - voidf(iel)) 217 vaporis = -cvap*min(0.d0, pressure(iel) - presat) & 218 *voidf(iel)*(1.d0 - voidf(iel)) 219 gamcav(iel) = condens + vaporis 220 if (gamcav(iel).lt.0) then 221 dgdpca(iel) = -cond*voidf(iel)*(1.d0 - voidf(iel)) 222 else 223 dgdpca(iel) = -cvap*voidf(iel)*(1.d0 - voidf(iel)) 224 endif 225 enddo 226 227 endif 228 229 end subroutine cavitation_compute_source_term 230 231 !============================================================================= 232 233 !> \brief Modify eddy viscosity using the Reboud correction: 234 !>\f[ 235 !> \mu_t'= \dfrac{\rho_v + (1-\alpha)^{mcav}(\rho_l-\rho_v)}{\rho}\mu_t. 236 !>\f] 237 !> 238 239 !> \param[in] crom density array 240 !> \param[in] voidf void fraction array 241 !> \param[in,out] visct turbulent viscosity 242 243 subroutine cavitation_correct_visc_turb (crom, voidf, visct) 244 245 use mesh 246 use numvar 247 use field 248 use vof 249 250 ! Arguments 251 252 double precision crom(ncelet), voidf(ncelet) 253 double precision visct(ncelet) 254 255 ! Local variables 256 257 integer iel 258 double precision frho 259 260 do iel = 1, ncel 261 frho = ( rho2 + (1.d0-voidf(iel))**mcav*(rho1 - rho2) ) & 262 /max(crom(iel),1.d-12) 263 visct(iel) = frho*visct(iel) 264 enddo 265 266 end subroutine cavitation_correct_visc_turb 267 268 !============================================================================= 269 270end module cavitation 271