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