From dc6d3756d8a9436a07a6951bc2a6995b4505f8d8 Mon Sep 17 00:00:00 2001 From: hepap <47506601+hepap@users.noreply.github.com> Date: Sat, 14 Oct 2023 12:19:42 -0400 Subject: [PATCH] Surface tension model for phase change (#904) Description of the problem Surface tension force was applied everywhere when treating solid-liquid phase change. This caused strange behaviour at the solid-liquid interface Description of the solution a surface tension model was added to compute the surface tension coefficient as: sigma=alpha_l*(sigma_0+dsigmadt*(T-T_0)), where alpha_l is the liquid fraction How Has This Been Tested? 2 tests added --- .../liquid_thermocapillary_droplet.output | 62 +++++ .../liquid_thermocapillary_droplet.prm | 237 ++++++++++++++++++ .../solid_thermocapillary_droplet.output | 62 +++++ .../solid_thermocapillary_droplet.prm | 237 ++++++++++++++++++ .../parameters/cfd/physical_properties.rst | 40 ++- include/core/parameters.h | 9 +- include/core/phase_change.h | 20 +- include/core/surface_tension_model.h | 161 +++++++++++- include/solvers/simulation_parameters.h | 34 ++- source/core/parameters.cc | 37 ++- source/core/surface_tension_model.cc | 4 + source/core/utilities.cc | 2 +- 12 files changed, 885 insertions(+), 20 deletions(-) create mode 100644 applications_tests/lethe-fluid/liquid_thermocapillary_droplet.output create mode 100644 applications_tests/lethe-fluid/liquid_thermocapillary_droplet.prm create mode 100644 applications_tests/lethe-fluid/solid_thermocapillary_droplet.output create mode 100644 applications_tests/lethe-fluid/solid_thermocapillary_droplet.prm diff --git a/applications_tests/lethe-fluid/liquid_thermocapillary_droplet.output b/applications_tests/lethe-fluid/liquid_thermocapillary_droplet.output new file mode 100644 index 0000000000..53cce6ca16 --- /dev/null +++ b/applications_tests/lethe-fluid/liquid_thermocapillary_droplet.output @@ -0,0 +1,62 @@ +Running on 1 MPI rank(s)... + Number of active cells: 4096 + Number of degrees of freedom: 12675 + Volume of triangulation: 36 + Number of thermal degrees of freedom: 4225 + Number of VOF degrees of freedom: 4225 + ++------------------------------------------+ +| VOF Barycenter | ++------------------------------------------+ + time x_vof y_vof vx_vof vy_vof +0.000000 -0.000000 -0.000000 0.000000 0.000000 + +******************************************************************************* +Transient iteration: 1 Time: 0.005 Time step: 0.005 CFL: 0 +******************************************************************************* + ++------------------------------------------+ +| VOF Barycenter | ++------------------------------------------+ + time x_vof y_vof vx_vof vy_vof +0.005000 -0.000000 -0.000000 0.001872 -0.000000 + +********************************************************************************** +Transient iteration: 2 Time: 0.01 Time step: 0.005 CFL: 0.000465832 +********************************************************************************** + ++------------------------------------------+ +| VOF Barycenter | ++------------------------------------------+ + time x_vof y_vof vx_vof vy_vof +0.010000 0.000009 -0.000000 0.003684 0.000000 + +********************************************************************************* +Transient iteration: 3 Time: 0.015 Time step: 0.005 CFL: 0.00084471 +********************************************************************************* + ++------------------------------------------+ +| VOF Barycenter | ++------------------------------------------+ + time x_vof y_vof vx_vof vy_vof +0.015000 0.000028 -0.000000 0.005439 0.000000 + +********************************************************************************* +Transient iteration: 4 Time: 0.02 Time step: 0.005 CFL: 0.00116533 +********************************************************************************* + ++------------------------------------------+ +| VOF Barycenter | ++------------------------------------------+ + time x_vof y_vof vx_vof vy_vof +0.020000 0.000055 -0.000000 0.007143 0.000000 + +********************************************************************************* +Transient iteration: 5 Time: 0.025 Time step: 0.005 CFL: 0.00144378 +********************************************************************************* + ++------------------------------------------+ +| VOF Barycenter | ++------------------------------------------+ + time x_vof y_vof vx_vof vy_vof +0.025000 0.000091 -0.000000 0.008798 0.000000 diff --git a/applications_tests/lethe-fluid/liquid_thermocapillary_droplet.prm b/applications_tests/lethe-fluid/liquid_thermocapillary_droplet.prm new file mode 100644 index 0000000000..7686b31013 --- /dev/null +++ b/applications_tests/lethe-fluid/liquid_thermocapillary_droplet.prm @@ -0,0 +1,237 @@ +#---------------------- +# Listing of Parameters +#---------------------- + +set dimension = 2 + +# Simulation and IO Control +#--------------------------------------------------- + +subsection simulation control + set method = bdf1 + set time end = 0.025 + set time step = 0.005 + set adapt = false + set max cfl = 0 + set output frequency = 0 +end + +#--------------------------------------------------- +# Multiphysics +#--------------------------------------------------- + +subsection multiphysics + set VOF = true + set heat transfer = true +end + +#--------------------------------------------------- +# VOF +#--------------------------------------------------- + +subsection VOF + subsection surface tension force + set enable = true + set phase fraction gradient diffusion factor = 4 + set curvature diffusion factor = 1 + set output auxiliary fields = true + set enable marangoni effect = true + end +end + +#--------------------------------------------------- +# Initial condition +#--------------------------------------------------- + +subsection initial conditions + set type = nodal + subsection uvwp + set Function expression = 0; 0; 0 + end + subsection VOF + set Function expression = if (x*x + y*y < 1.0*1.0 , 1, 0) + subsection projection step + set enable = true + set diffusion factor = 1.0 + end + end + subsection temperature + set Function expression = x+3 + end +end + +#--------------------------------------------------- +# Physical Properties +#--------------------------------------------------- + +subsection physical properties + set number of fluids = 2 + subsection fluid 1 + set density = 1 + set kinematic viscosity = 1 + set thermal conductivity = 100 + end + subsection fluid 0 + set density = 1 + set kinematic viscosity = 1 + set thermal conductivity = 100 + end + + set number of material interactions = 1 + subsection material interaction 0 + set type = fluid-fluid + subsection fluid-fluid interaction + set first fluid id = 0 + set second fluid id = 1 + set surface tension model = phase change + set surface tension coefficient = 3.0 + set reference state temperature = 0.0 + set temperature-driven surface tension gradient = -1.0 + set liquidus temperature = 0 + set solidus temperature = -1 + end + end +end + +#--------------------------------------------------- +# Mesh +#--------------------------------------------------- + +subsection mesh + set type = dealii + set grid type = subdivided_hyper_rectangle + set grid arguments = 1, 1: -3.0, -3.0: 3.0, 3.0: true + set initial refinement = 6 +end + +# -------------------------------------------------- +# Boundary Conditions +#--------------------------------------------------- + +subsection boundary conditions + set number = 4 + subsection bc 0 + set id = 0 + set type = noslip + end + subsection bc 1 + set id = 1 + set type = noslip + end + subsection bc 2 + set id = 2 + set type = noslip + end + subsection bc 3 + set id = 3 + set type = noslip + end +end + +subsection boundary conditions heat transfer + set number = 4 + subsection bc 0 + set id = 0 + set type = temperature + set value = 0 + end + subsection bc 1 + set id = 1 + set type = temperature + set value = 6 + end + subsection bc 2 + set id = 2 + set type = noflux + end + subsection bc 3 + set id = 3 + set type = noflux + end +end + +#--------------------------------------------------- +# Post-processing +#--------------------------------------------------- + +subsection post-processing + set verbosity = verbose + set calculate barycenter = true +end + + +#--------------------------------------------------- +# FEM +#--------------------------------------------------- + +subsection FEM + set velocity order = 1 + set pressure order = 1 + set VOF order = 1 + set temperature order = 1 +end + +# -------------------------------------------------- +# Non-Linear Solver Control +#--------------------------------------------------- + +subsection non-linear solver + subsection VOF + set tolerance = 1e-5 + set max iterations = 100 + set verbosity = quiet + end + subsection heat transfer + set tolerance = 1e-5 + set max iterations = 100 + set verbosity = quiet + end + subsection fluid dynamics + set tolerance = 1e-5 + set max iterations = 100 + set verbosity = quiet + end +end + +#--------------------------------------------------- +# Linear Solver Control +#--------------------------------------------------- + +subsection linear solver + subsection fluid dynamics + set verbosity = quiet + set method = gmres + set max iters = 8000 + set relative residual = 1e-3 + set minimum residual = 1e-7 + set preconditioner = ilu + set ilu preconditioner fill = 0 + set ilu preconditioner absolute tolerance = 1e-12 + set ilu preconditioner relative tolerance = 1.00 + set max krylov vectors = 200 + end + subsection VOF + set verbosity = quiet + set method = gmres + set max iters = 8000 + set relative residual = 1e-3 + set minimum residual = 1e-9 + set preconditioner = ilu + set ilu preconditioner fill = 0 + set ilu preconditioner absolute tolerance = 1e-12 + set ilu preconditioner relative tolerance = 1.00 + set max krylov vectors = 200 + end + subsection heat transfer + set verbosity = quiet + set method = gmres + set max iters = 8000 + set relative residual = 1e-3 + set minimum residual = 1e-9 + set preconditioner = ilu + set ilu preconditioner fill = 0 + set ilu preconditioner absolute tolerance = 1e-12 + set ilu preconditioner relative tolerance = 1.00 + set max krylov vectors = 200 + end +end diff --git a/applications_tests/lethe-fluid/solid_thermocapillary_droplet.output b/applications_tests/lethe-fluid/solid_thermocapillary_droplet.output new file mode 100644 index 0000000000..b95420a6cb --- /dev/null +++ b/applications_tests/lethe-fluid/solid_thermocapillary_droplet.output @@ -0,0 +1,62 @@ +Running on 1 MPI rank(s)... + Number of active cells: 4096 + Number of degrees of freedom: 12675 + Volume of triangulation: 36 + Number of thermal degrees of freedom: 4225 + Number of VOF degrees of freedom: 4225 + ++------------------------------------------+ +| VOF Barycenter | ++------------------------------------------+ + time x_vof y_vof vx_vof vy_vof +0.000000 -0.000000 -0.000000 0.000000 0.000000 + +******************************************************************************* +Transient iteration: 1 Time: 0.005 Time step: 0.005 CFL: 0 +******************************************************************************* + ++------------------------------------------+ +| VOF Barycenter | ++------------------------------------------+ + time x_vof y_vof vx_vof vy_vof +0.005000 -0.000000 -0.000000 0.000000 0.000000 + +******************************************************************************* +Transient iteration: 2 Time: 0.01 Time step: 0.005 CFL: 0 +******************************************************************************* + ++------------------------------------------+ +| VOF Barycenter | ++------------------------------------------+ + time x_vof y_vof vx_vof vy_vof +0.010000 -0.000000 -0.000000 0.000000 0.000000 + +******************************************************************************* +Transient iteration: 3 Time: 0.015 Time step: 0.005 CFL: 0 +******************************************************************************* + ++------------------------------------------+ +| VOF Barycenter | ++------------------------------------------+ + time x_vof y_vof vx_vof vy_vof +0.015000 -0.000000 -0.000000 0.000000 0.000000 + +******************************************************************************* +Transient iteration: 4 Time: 0.02 Time step: 0.005 CFL: 0 +******************************************************************************* + ++------------------------------------------+ +| VOF Barycenter | ++------------------------------------------+ + time x_vof y_vof vx_vof vy_vof +0.020000 -0.000000 -0.000000 0.000000 0.000000 + +******************************************************************************* +Transient iteration: 5 Time: 0.025 Time step: 0.005 CFL: 0 +******************************************************************************* + ++------------------------------------------+ +| VOF Barycenter | ++------------------------------------------+ + time x_vof y_vof vx_vof vy_vof +0.025000 -0.000000 -0.000000 0.000000 0.000000 diff --git a/applications_tests/lethe-fluid/solid_thermocapillary_droplet.prm b/applications_tests/lethe-fluid/solid_thermocapillary_droplet.prm new file mode 100644 index 0000000000..3ea53bd234 --- /dev/null +++ b/applications_tests/lethe-fluid/solid_thermocapillary_droplet.prm @@ -0,0 +1,237 @@ +#---------------------- +# Listing of Parameters +#---------------------- + +set dimension = 2 + +# Simulation and IO Control +#--------------------------------------------------- + +subsection simulation control + set method = bdf1 + set time end = 0.025 + set time step = 0.005 + set adapt = false + set max cfl = 0 + set output frequency = 0 +end + +#--------------------------------------------------- +# Multiphysics +#--------------------------------------------------- + +subsection multiphysics + set VOF = true + set heat transfer = true +end + +#--------------------------------------------------- +# VOF +#--------------------------------------------------- + +subsection VOF + subsection surface tension force + set enable = true + set phase fraction gradient diffusion factor = 4 + set curvature diffusion factor = 1 + set output auxiliary fields = true + set enable marangoni effect = true + end +end + +#--------------------------------------------------- +# Initial condition +#--------------------------------------------------- + +subsection initial conditions + set type = nodal + subsection uvwp + set Function expression = 0; 0; 0 + end + subsection VOF + set Function expression = if (x*x + y*y < 1.0*1.0 , 1, 0) + subsection projection step + set enable = true + set diffusion factor = 1.0 + end + end + subsection temperature + set Function expression = x+3 + end +end + +#--------------------------------------------------- +# Physical Properties +#--------------------------------------------------- + +subsection physical properties + set number of fluids = 2 + subsection fluid 1 + set density = 1 + set kinematic viscosity = 1 + set thermal conductivity = 100 + end + subsection fluid 0 + set density = 1 + set kinematic viscosity = 1 + set thermal conductivity = 100 + end + + set number of material interactions = 1 + subsection material interaction 0 + set type = fluid-fluid + subsection fluid-fluid interaction + set first fluid id = 0 + set second fluid id = 1 + set surface tension model = phase change + set surface tension coefficient = 3.0 + set reference state temperature = 0.0 + set temperature-driven surface tension gradient = -1.0 + set liquidus temperature = 7 + set solidus temperature = 6 + end + end +end + +#--------------------------------------------------- +# Mesh +#--------------------------------------------------- + +subsection mesh + set type = dealii + set grid type = subdivided_hyper_rectangle + set grid arguments = 1, 1: -3.0, -3.0: 3.0, 3.0: true + set initial refinement = 6 +end + +# -------------------------------------------------- +# Boundary Conditions +#--------------------------------------------------- + +subsection boundary conditions + set number = 4 + subsection bc 0 + set id = 0 + set type = noslip + end + subsection bc 1 + set id = 1 + set type = noslip + end + subsection bc 2 + set id = 2 + set type = noslip + end + subsection bc 3 + set id = 3 + set type = noslip + end +end + +subsection boundary conditions heat transfer + set number = 4 + subsection bc 0 + set id = 0 + set type = temperature + set value = 0 + end + subsection bc 1 + set id = 1 + set type = temperature + set value = 6 + end + subsection bc 2 + set id = 2 + set type = noflux + end + subsection bc 3 + set id = 3 + set type = noflux + end +end + +#--------------------------------------------------- +# Post-processing +#--------------------------------------------------- + +subsection post-processing + set verbosity = verbose + set calculate barycenter = true +end + + +#--------------------------------------------------- +# FEM +#--------------------------------------------------- + +subsection FEM + set velocity order = 1 + set pressure order = 1 + set VOF order = 1 + set temperature order = 1 +end + +# -------------------------------------------------- +# Non-Linear Solver Control +#--------------------------------------------------- + +subsection non-linear solver + subsection VOF + set tolerance = 1e-5 + set max iterations = 100 + set verbosity = quiet + end + subsection heat transfer + set tolerance = 1e-5 + set max iterations = 100 + set verbosity = quiet + end + subsection fluid dynamics + set tolerance = 1e-5 + set max iterations = 100 + set verbosity = quiet + end +end + +#--------------------------------------------------- +# Linear Solver Control +#--------------------------------------------------- + +subsection linear solver + subsection fluid dynamics + set verbosity = quiet + set method = gmres + set max iters = 8000 + set relative residual = 1e-3 + set minimum residual = 1e-7 + set preconditioner = ilu + set ilu preconditioner fill = 0 + set ilu preconditioner absolute tolerance = 1e-12 + set ilu preconditioner relative tolerance = 1.00 + set max krylov vectors = 200 + end + subsection VOF + set verbosity = quiet + set method = gmres + set max iters = 8000 + set relative residual = 1e-3 + set minimum residual = 1e-9 + set preconditioner = ilu + set ilu preconditioner fill = 0 + set ilu preconditioner absolute tolerance = 1e-12 + set ilu preconditioner relative tolerance = 1.00 + set max krylov vectors = 200 + end + subsection heat transfer + set verbosity = quiet + set method = gmres + set max iters = 8000 + set relative residual = 1e-3 + set minimum residual = 1e-9 + set preconditioner = ilu + set ilu preconditioner fill = 0 + set ilu preconditioner absolute tolerance = 1e-12 + set ilu preconditioner relative tolerance = 1.00 + set max krylov vectors = 200 + end +end diff --git a/doc/source/parameters/cfd/physical_properties.rst b/doc/source/parameters/cfd/physical_properties.rst index 482376a122..8593d7cedc 100644 --- a/doc/source/parameters/cfd/physical_properties.rst +++ b/doc/source/parameters/cfd/physical_properties.rst @@ -45,8 +45,12 @@ Physical Properties set second fluid id = 1 # Surface tension - set surface tension model = constant - set surface tension coefficient = 0 + set surface tension model = constant + set surface tension coefficient = 0 + set reference state temperature = 0 + set temperature-driven surface tension gradient = 0 + set liquidus temperature = 0 + set solidus temperature = 0 # Mobility Cahn-Hilliard set cahn hilliard mobility model = constant @@ -63,6 +67,8 @@ Physical Properties set surface tension coefficient = 0 set reference state temperature = 0 set temperature-driven surface tension gradient = 0 + set liquidus temperature = 0 + set solidus temperature = 0 end end end @@ -114,7 +120,7 @@ Physical Properties .. attention:: The ``second fluid id`` should be greater than the ``first fluid id``. - * The ``surface tension model`` specifies the model used to calculate the surface tension coefficient of the fluid-fluid pair. At the moment, a ``constant`` and a ``linear`` model are supported. For more detail on the surface tension models, see `Surface Tension Models`_. + * The ``surface tension model`` specifies the model used to calculate the surface tension coefficient of the fluid-fluid pair. At the moment, ``constant``, ``linear``, and ``phase change`` models are supported. For more details on the surface tension models, see `Surface Tension Models`_. * The ``surface tension coefficient`` parameter is a constant surface tension coefficient of the two interacting fluids in units of :math:`\text{Mass} \cdot \text{Time}^{-2}`. In SI, this is :math:`\text{N} \cdot \text{m}^{-1}`. The surface tension coefficient is used as defined in the Weber number (:math:`We`): @@ -126,6 +132,8 @@ Physical Properties * The ``reference state temperature`` parameter is the temperature of the reference state at which the ``surface tension coefficient`` is evaluated. This parameter is used in the calculation of the surface tension using the ``linear`` surface tension model (see `Surface Tension Models`_). * The ``temperature-driven surface tension gradient`` parameter is the surface tension gradient with respect to the temperature of the two interacting fluids in units of :math:`\text{Mass} \cdot \text{Time}^{-2} \cdot \text{Temperature}^{-1}`. In SI, this is :math:`\text{N} \cdot \text{m}^{-1} \cdot \text{K}^{-1}`. This parameter is used in the calculation of the surface tension using the ``linear`` surface tension model (see `Surface Tension Models`_). + + * The ``solidus temperature`` and ``liquidus temperature`` parameters are used in the calculation of the surface tension using the ``phase change`` surface tension model (see `Surface Tension Models`_). * The ``cahn hilliard mobility model`` specifies the model used to calculate the mobility used in the Cahn-Hilliard equations for the pair of fluid. Two models are available: a ``constant`` mobility and a ``quartic`` mobility. The reader is refered to :doc:`cahn_hilliard` for more details. @@ -500,15 +508,35 @@ Interface Physical Property Models Surface Tension Models ~~~~~~~~~~~~~~~~~~~~~~~ -Lethe supports two types of surface tension models: ``constant`` and ``linear``. A ``constant`` surface tension model assumes a constant value of surface tension, while a ``linear`` surface tension assumes that the surface tension evolves linearly with the temperature: +Lethe supports three types of surface tension models: ``constant``, ``linear``, and ``phase change``. A ``constant`` surface tension model assumes a constant value of surface tension, while a ``linear`` surface tension assumes that the surface tension evolves linearly with the temperature: .. math:: \sigma(T) = \sigma_0 + \frac{d\sigma}{dT} (T-T_0) where :math:`\sigma_0` is the ``surface tension coefficient`` evaluated at ``reference state temperature`` :math:`T_0` and :math:`\frac{d\sigma}{dT}` is the ``surface tension gradient`` with respect to the temperature :math:`T`. +For problems treating solid-liquid phase change, the ``phase change`` model is intended to apply the surface tension force only when the fluid is liquid such that: + +.. math:: + \sigma(T) = + \begin{cases} + 0 &\quad\text{if}\; T`_ B. Blais and F. Ilinca, “Development and validation of a stabilized immersed boundary CFD model for freezing and melting with natural convection,” *Comput. Fluids*, vol. 172, pp. 564–581, Aug. 2018, doi: 10.1016/j.compfluid.2018.03.037. `[2] `_ E. Bretin, R. Denis, S. Masnou, A. Sengers, and G. Terii, “A multiphase Cahn-Hilliard system with mobilities and the numerical simulation of dewetting.” arXiv, Apr. 18, 2023. doi: 10.48550/arXiv.2105.09627. - - diff --git a/include/core/parameters.h b/include/core/parameters.h index e9e895f076..d8ccee279d 100644 --- a/include/core/parameters.h +++ b/include/core/parameters.h @@ -304,6 +304,12 @@ namespace Parameters // N/(m*K) double surface_tension_gradient; + // Solidus temperature - Units in K + double T_solidus; + + // Liquidus temperature - Units in K + double T_liquidus; + void declare_parameters(ParameterHandler &prm); void @@ -423,7 +429,8 @@ namespace Parameters enum class SurfaceTensionModel { constant, - linear + linear, + phase_change } surface_tension_model; SurfaceTensionParameters surface_tension_parameters; diff --git a/include/core/phase_change.h b/include/core/phase_change.h index 6727056a31..12574f18d7 100644 --- a/include/core/phase_change.h +++ b/include/core/phase_change.h @@ -25,7 +25,7 @@ using namespace dealii; * @brief calculate_liquid_fraction Calculates the liquid fraction of a phase change material at * a temperature T * - * @param T temperature at which to calculate the solid fraction + * @param T temperature at which to calculate the liquid fraction * @return value of the liquid_fraction * */ @@ -41,4 +41,22 @@ calculate_liquid_fraction( 1.); } +/** + * @brief calculate_liquid_fraction Calculates the liquid fraction of a phase change material at + * a temperature T + * + * @param T temperature at which to calculate the liquid fraction + * @param solidus temperature + * @param liquidus temperature + * @return value of the liquid_fraction + * + */ +inline double +calculate_liquid_fraction(const double &T, + const double T_solidus, + const double T_liquidus) +{ + return std::min(std::max((T - T_solidus) / (T_liquidus - T_solidus), 0.), 1.); +} + #endif diff --git a/include/core/surface_tension_model.h b/include/core/surface_tension_model.h index 531c3a5f00..cd1e6a89e3 100644 --- a/include/core/surface_tension_model.h +++ b/include/core/surface_tension_model.h @@ -18,6 +18,7 @@ #define lethe_surface_tension_model_h #include +#include /** * @brief SurfaceTensionModel. Abstract class that allows to calculate the @@ -145,6 +146,9 @@ class SurfaceTensionLinear : public SurfaceTensionModel public: /** * @brief Default constructor + * @param p_surface_tension_parameters The parameters for the surface tension + * model, including the surface tension coefficient sigma_0 and surface + * tension gradient dsimga/dT at the reference temperature T_0 */ SurfaceTensionLinear( const Parameters::SurfaceTensionParameters &p_surface_tension_parameters) @@ -231,10 +235,165 @@ class SurfaceTensionLinear : public SurfaceTensionModel std::fill(jacobian_vector.begin(), jacobian_vector.end(), 0); } -private: +protected: const double surface_tension_coefficient; const double T_0; const double surface_tension_gradient; }; +/** + * @brief Linear surface tension with phase change. The surface tension is given by: + * alpha_l*(sigma_0 + dsigma/dT * (T-T_0)). + */ +class SurfaceTensionPhaseChange : public SurfaceTensionLinear +{ +public: + /** + * @brief Default constructor + * @param p_surface_tension_parameters The parameters for the surface tension + * model, including the surface tension coefficient sigma_0 and surface + * tension gradient dsimga/dT at the reference temperature T_0, and the + * solidus and liquidus temperatures T_solidus and T_liquidus respectively. + */ + SurfaceTensionPhaseChange( + const Parameters::SurfaceTensionParameters &p_surface_tension_parameters) + : SurfaceTensionLinear(p_surface_tension_parameters) + , T_solidus(p_surface_tension_parameters.T_solidus) + , T_liquidus(p_surface_tension_parameters.T_liquidus) + { + this->model_depends_on[field::temperature] = true; + } + + /** + * @brief value Calculates the surface tension coefficient. + * @param field_vectors Vectors of the fields on which the surface tension + * depends. In this case, it depends on the temperature. + * @return value of the physical property calculated with the fields_value. + */ + double + value(const std::map &fields_value) override + { + const double temperature = fields_value.at(field::temperature); + + double surface_tension; + + if (temperature < T_solidus) + surface_tension = 0.0; + else if (temperature > T_liquidus) + surface_tension = surface_tension_coefficient + + surface_tension_gradient * (temperature - T_0); + else + { + const double l_frac = + calculate_liquid_fraction(temperature, T_solidus, T_liquidus); + + surface_tension = (surface_tension_coefficient + + surface_tension_gradient * (temperature - T_0)) * + l_frac; + } + return surface_tension; + } + + /** + * @brief vector_value Calculates the vector of surface tension coefficient. + * @param field_vectors Vectors of the fields on which the surface tension + * depends. In this case, it depends on the temperature. + * @param property_vector Vectors of the surface tension coefficient values. + */ + void + vector_value(const std::map> &field_vectors, + std::vector &property_vector) override + { + const std::vector &temperature = + field_vectors.at(field::temperature); + for (unsigned int i = 0; i < property_vector.size(); ++i) + if (temperature[i] < T_solidus) + property_vector[i] = 0.0; + else if (temperature[i] > T_liquidus) + property_vector[i] = surface_tension_coefficient + + surface_tension_gradient * (temperature[i] - T_0); + else + { + const double l_frac = + calculate_liquid_fraction(temperature[i], T_solidus, T_liquidus); + + property_vector[i] = + (surface_tension_coefficient + + surface_tension_gradient * (temperature[i] - T_0)) * + l_frac; + } + } + + /** + * @brief jacobian Calculates the jacobian (the partial derivative) of the + * surface tension coefficient with respect to a field. + * @param field_values Value of the various fields on which the property may + * depend. + * @param id Indicator of the field with respect to which the jacobian + * should be calculated. + * @return value of the partial derivative of the surface tension coefficient + * with respect to the field. + */ + double + jacobian(const std::map &field_values, field id) override + { + if (id == field::temperature) + { + const double temperature = field_values.at(field::temperature); + if (temperature < T_solidus) + return 0; + else if (temperature > T_liquidus) + return surface_tension_gradient; + else + { + const double l_frac = + calculate_liquid_fraction(temperature, T_solidus, T_liquidus); + return surface_tension_gradient * l_frac; + } + } + else + return 0; + } + + /** + * @brief vector_jacobian Calculates the derivative of the surface tension + * coefficient with respect to a field. + * @param field_vectors Vector for the values of the fields used to evaluate + * the property. + * @param id Identifier of the field with respect to which a derivative should + * be calculated. + * @param jacobian vector of the value of the derivative of the surface + * tension coefficient with respect to the field id. + */ + void + vector_jacobian(const std::map> &field_vectors, + const field id, + std::vector &jacobian_vector) override + { + if (id == field::temperature) + { + const std::vector &temperature = + field_vectors.at(field::temperature); + for (unsigned int i = 0; i < jacobian_vector.size(); ++i) + if (temperature[i] < T_solidus) + jacobian_vector[i] = 0.0; + else if (temperature[i] > T_liquidus) + jacobian_vector[i] = surface_tension_gradient; + else + { + const double l_frac = calculate_liquid_fraction(temperature[i], + T_solidus, + T_liquidus); + jacobian_vector[i] = surface_tension_gradient * l_frac; + } + } + else + std::fill(jacobian_vector.begin(), jacobian_vector.end(), 0); + } + +private: + const double T_solidus; + const double T_liquidus; +}; + #endif diff --git a/include/solvers/simulation_parameters.h b/include/solvers/simulation_parameters.h index 95c26c7822..8f93b13c30 100644 --- a/include/solvers/simulation_parameters.h +++ b/include/solvers/simulation_parameters.h @@ -221,6 +221,16 @@ class SimulationParameters " set temperature-driven surface tension gradient = $value_of_gradient\n" " end\n"); + std::string phase_change_surface_tension_model( + " subsection fluid-fluid interaction\n" + " set first fluid id = 0\n" + " set second fluid id = 1\n" + " set surface tension model = phase change\n" + " set surface tension coefficient = $value_of_coefficient\n" + " set temperature-driven surface tension gradient = $value_of_gradient\n" + " set solidus temperature = $value_of_solidus_temperature\n" + " set liquidus temperature = $value_of_liquidus_temperature\n" + " end\n"); if (!multiphysics.vof_parameters.surface_tension_force .enable_marangoni_effect) // constant surface tension model { @@ -287,7 +297,13 @@ class SimulationParameters " set number of material interactions = 1\n" " subsection material interaction 0\n" " set type = fluid-fluid\n" + - linear_surface_tension_model + " end\n"); + linear_surface_tension_model + + "\n" + "or: \n\n" + " set number of material interactions = 1\n" + " subsection material interaction 0\n" + " set type = fluid-fluid\n" + + phase_change_surface_tension_model + " end\n"); } else { @@ -302,7 +318,13 @@ class SimulationParameters "effect. In subsection physical properties, use:\n\n" " subsection material interaction $material_interaction_id\n" " set type = fluid-fluid\n" + - linear_surface_tension_model + " end\n"); + linear_surface_tension_model + + "\n" + "or:\n\n" + " set number of material interactions = 1\n" + " subsection material interaction 0\n" + " set type = fluid-fluid\n" + + phase_change_surface_tension_model + " end\n"); } if (is_constant_surface_tension_model( physical_properties.material_interactions)) @@ -315,7 +337,13 @@ class SimulationParameters "effect. In subsection physical properties, use:\n\n" " subsection material interaction $material_interaction_id\n" " set type = fluid-fluid\n" + - linear_surface_tension_model + " end\n"); + linear_surface_tension_model + + "\n" + "or:\n\n" + " set number of material interactions = 1\n" + " subsection material interaction 0\n" + " set type = fluid-fluid\n" + + phase_change_surface_tension_model + " end\n"); } } } diff --git a/source/core/parameters.cc b/source/core/parameters.cc index efc8065cdc..e4149a773d 100644 --- a/source/core/parameters.cc +++ b/source/core/parameters.cc @@ -431,6 +431,16 @@ namespace Parameters "0.0", Patterns::Double(), "Surface tension gradient with respect to the temperature for the corresponding pair of fluids or fluid-solid pair"); + prm.declare_entry( + "solidus temperature", + "0", + Patterns::Double(), + "Temperature of the solidus for the corresponding pair of fluids or fluid-solid pair"); + prm.declare_entry( + "liquidus temperature", + "1", + Patterns::Double(), + "Temperature of the liquidus for the corresponding pair of fluids or fluid-solid pair"); } void @@ -440,6 +450,11 @@ namespace Parameters T_0 = prm.get_double("reference state temperature"); surface_tension_gradient = prm.get_double("temperature-driven surface tension gradient"); + T_solidus = prm.get_double("solidus temperature"); + T_liquidus = prm.get_double("liquidus temperature"); + + Assert(T_liquidus > T_solidus, + PhaseChangeIntervalError(T_liquidus, T_solidus)); } void @@ -985,9 +1000,9 @@ namespace Parameters prm.declare_entry( "surface tension model", "constant", - Patterns::Selection("constant|linear"), + Patterns::Selection("constant|linear|phase change"), "Model used for the calculation of the surface tension coefficient\n" - "The choices are ."); + "The choices are ."); surface_tension_parameters.declare_parameters(prm); // Cahn-Hilliard mobility @@ -1017,9 +1032,9 @@ namespace Parameters prm.declare_entry( "surface tension model", "constant", - Patterns::Selection("constant|linear"), + Patterns::Selection("constant|linear|phase change"), "Model used for the calculation of the surface tension coefficient\n" - "The choices are ."); + "The choices are ."); surface_tension_parameters.declare_parameters(prm); } prm.leave_subsection(); @@ -1069,9 +1084,14 @@ namespace Parameters surface_tension_model = SurfaceTensionModel::linear; surface_tension_parameters.parse_parameters(prm); } + else if (op == "phase change") + { + surface_tension_model = SurfaceTensionModel::phase_change; + surface_tension_parameters.parse_parameters(prm); + } else throw(std::runtime_error( - "Invalid surface tension model. The choices are .")); + "Invalid surface tension model. The choices are .")); // Cahn-Hilliard mobility op = prm.get("cahn hilliard mobility model"); @@ -1114,9 +1134,14 @@ namespace Parameters surface_tension_model = SurfaceTensionModel::linear; surface_tension_parameters.parse_parameters(prm); } + else if (op == "phase change") + { + surface_tension_model = SurfaceTensionModel::phase_change; + surface_tension_parameters.parse_parameters(prm); + } else throw(std::runtime_error( - "Invalid surface tension model. The choices are .")); + "Invalid surface tension model. The choices are .")); std::pair, SurfaceTensionModel> fluid_solid_surface_tension_interaction(fluid_solid_interaction, surface_tension_model); diff --git a/source/core/surface_tension_model.cc b/source/core/surface_tension_model.cc index 7e4f20aedb..220129611e 100644 --- a/source/core/surface_tension_model.cc +++ b/source/core/surface_tension_model.cc @@ -24,6 +24,10 @@ SurfaceTensionModel::model_cast( Parameters::MaterialInteractions::SurfaceTensionModel::linear) return std::make_shared( material_interaction_parameters.surface_tension_parameters); + else if (material_interaction_parameters.surface_tension_model == + Parameters::MaterialInteractions::SurfaceTensionModel::phase_change) + return std::make_shared( + material_interaction_parameters.surface_tension_parameters); else return std::make_shared( material_interaction_parameters.surface_tension_parameters diff --git a/source/core/utilities.cc b/source/core/utilities.cc index 26dca6714f..c8a9892b60 100644 --- a/source/core/utilities.cc +++ b/source/core/utilities.cc @@ -43,7 +43,7 @@ make_table_scalars_tensors( const std::vector &dependent_column_name, const unsigned int display_precision) { - AssertDimension(dependent_column_name.size(), 3 * dim); + AssertDimension(dependent_column_name.size(), dependent_vectors.size() * dim); TableHandler table; unsigned int vect_index = 0;