diff --git a/examples/reactmicp/unsaturated/acc_carbo.cpp b/examples/reactmicp/unsaturated/acc_carbo.cpp index e67f9f0..603b600 100644 --- a/examples/reactmicp/unsaturated/acc_carbo.cpp +++ b/examples/reactmicp/unsaturated/acc_carbo.cpp @@ -1,521 +1,508 @@ /*------------------------------------------------------------------------------ Copyright (c)2015 F. Georget , Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ------------------------------------------------------------------------------*/ //! \file acc_carbo.cpp //! \brief Accelerated carbonation example //! //! This is an example to show the features of the unsaturated //! reactive transport solver //! It is a simple study of the accelerated carbonation of cement paste //! (with drying) //! //! #include "reactmicp_unsaturated.hpp" #include "reactmicp/systems/unsaturated/variables.hpp" #include #include "reactmicp/io/hdf5_unsaturated.hpp" #include "utils/io/specmicp_hdf5.hpp" #include "dfpm/io/hdf5_mesh.hpp" #include "database/io/hdf5_database.hpp" #include "dfpmsolver/parabolic_structs.hpp" +#include "physics/unsaturated_laws.hpp" + #include using namespace specmicp; using namespace specmicp::reactmicp; using namespace specmicp::reactmicp::systems::unsaturated; using VariablesBase = solver::VariablesBase; using VariablesBasePtr = solver::VariablesBasePtr; using StaggerReturnCode = solver::StaggerReturnCode; int main(); class UpscalingStagger; database::RawDatabasePtr get_database(); AdimensionalSystemSolution get_initial_condition( database::RawDatabasePtr the_database, const units::UnitsSet& units_set ); // Upscaling stagger -> user models class UpscalingStagger: public solver::UpscalingStaggerBase { public: UpscalingStagger() {} ~UpscalingStagger() {} //! \brief Initialize the stagger at the beginning of the computation //! //! \param var a shared_ptr to the variables void initialize(VariablesBasePtr var); //! \brief Initialize the stagger at the beginning of an iteration //! //! This is where the predictor can be saved, the first trivial iteration done, ... //! //! \param dt the duration of the timestep //! \param var a shared_ptr to the variables void initialize_timestep(scalar_t dt, VariablesBasePtr var); //! \brief Solve the equation for the timestep //! //! \param var a shared_ptr to the variables StaggerReturnCode restart_timestep(VariablesBasePtr var); // User models scalar_t capillary_pressure(index_t node, scalar_t saturation); scalar_t vapor_pressure(index_t node, scalar_t saturation); scalar_t relative_liquid_diffusivity(index_t node, scalar_t saturation); scalar_t relative_liquid_permeability(index_t node, scalar_t saturation); scalar_t relative_gas_diffusivity(index_t node, scalar_t saturation); - scalar_t ln_hr(scalar_t saturation); - scalar_t intrinsic_permeability(scalar_t porosity); scalar_t intrinsic_liquid_diffusivity(scalar_t porosity); scalar_t resistance_gas_diffusivity(scalar_t porosity); private: + scalar_t m_a {37.5479e6}; + scalar_t m_b {2.1684}; scalar_t m_alpha {7.27}; scalar_t m_beta {0.440}; scalar_t m_pv_sat {3170}; scalar_t m_exp_r_l_d {3.00}; scalar_t m_exp_i_g_d {2.74}; scalar_t m_exp_r_g_d {4.20}; scalar_t m_exp_r_l_k {0.4191}; scalar_t m_k_0 {1e-21*1e4}; scalar_t m_d_0 {2.3e-13*1e4}; scalar_t m_phi_0 {0.2}; }; // implementations int main() { init_logger(&std::cout, logger::Warning); index_t nb_nodes = 50; scalar_t dx = 0.05; // cm scalar_t cross_section = 5; // cm^2 units::UnitsSet units_set; units_set.length = units::LengthUnit::centimeter; specmicp::RawDatabasePtr the_database = get_database(); AdimensionalSystemSolution init = get_initial_condition( the_database, units_set); AdimensionalSystemSolutionExtractor extr(init, the_database, units_set); std::cout << extr.porosity() << " - " << extr.saturation_water() << std::endl; //return EXIT_SUCCESS; mesh::Mesh1DPtr the_mesh = mesh::uniform_mesh1d({dx, nb_nodes, cross_section}); index_t id_co2 = the_database->get_id_component("CO2"); VariablesInterface variable_interface(the_mesh, the_database, {0, id_co2}, units_set); variable_interface.set_binary_gas_diffusivity(0, 0.240); variable_interface.set_binary_gas_diffusivity(id_co2, 0.160); std::shared_ptr upscaling_stagger = std::make_shared(); auto* call_ptr = upscaling_stagger.get(); variable_interface.set_vapor_pressure_model( [call_ptr](index_t node, scalar_t sat) -> scalar_t { return call_ptr->vapor_pressure(node, sat); }); for (index_t node=1; nodenb_nodes(); ++node) { variable_interface.initialize_variables(node, extr); } variable_interface.set_porosity(0, 1.0); variable_interface.set_liquid_saturation(0, 0.0); variable_interface.set_partial_pressure(id_co2, 0, 100.0); variable_interface.set_partial_pressure(0, 0, 1000.0); std::shared_ptr vars = variable_interface.get_variables(); vars->set_aqueous_scaling(0, 1.0e3/18.3e-3*1e-6); for (auto aqc: the_database->range_aqueous_component()) { vars->set_aqueous_scaling(aqc, 100e-6); } //vars->set_gaseous_scaling(0, 1.0/vars->get_rt()); //vars->set_gaseous_scaling(id_co2, 1.0/vars->get_rt()); std::cout << "Youpla : " << vars->get_rt() << std::endl; std::shared_ptr eq_constraints = equilibrium_constraints(the_mesh); eq_constraints->add_fixed_node(0); eq_constraints->set_constraint_for_all(eq_constraints->new_constraints()); eq_constraints->get_options().units_set = units_set; micpsolver::MiCPSolverOptions& copts = eq_constraints->get_options().solver_options; copts.set_maximum_iterations(200); copts.set_maximum_step_length(10, 200); copts.set_tolerance(1e-10, 1e-16); eq_constraints->get_options().system_options.non_ideality_tolerance = 1e-14; std::shared_ptr chemistry_stagger = equilibrium_stagger(vars, eq_constraints); upscaling_stagger->initialize(vars); vars->set_relative_variables(0); chemistry_stagger->initialize(vars); TransportConstraints transport_constraints; transport_constraints.add_gas_node(0); transport_constraints.add_fixed_node(0); std::shared_ptr transport_stagger = unsaturated_transport_stagger(vars, transport_constraints, true, true); transport_stagger->initialize(vars); ((*eq_constraints)[0]).set_water_partial_pressure_model( std::bind(&UpscalingStagger::vapor_pressure, upscaling_stagger.get(), 1, std::placeholders::_1) ); auto check = variable_interface.check_variables(); if (check >= VariablesValidity::error) { throw std::runtime_error("Invalid variables"); } solver::ReactiveTransportSolver react_solver(transport_stagger, chemistry_stagger, upscaling_stagger); solver::ReactiveTransportOptions& opts = react_solver.get_options(); opts.residuals_tolerance = 1e-2; opts.good_enough_tolerance = 0.99; opts.maximum_iterations = 500; transport_stagger->get_saturation_options()->maximum_iterations = 5000; transport_stagger->get_saturation_options()->step_tolerance = 1e-16; transport_stagger->get_saturation_options()->residuals_tolerance = 1e-10; //transport_stagger->get_gas_options(3)->maximum_iterations = 5000; opts.step_tolerance = 1e-10; io::HDF5File save_file("out_acc_carbo.hdf5", io::HDF5_OpenMode::CreateTruncate); io::save_database_labels(save_file, "database", "/", the_database); io::save_mesh_coordinates(save_file, "mesh", "/", the_mesh); io::UnsaturatedHDF5Saver saver(vars); saver.save_timepoint(save_file, "0", "/"); save_file.close(); solver::SimulationInformation simul_info("acc_carbo", 900); auto out_pol = [&saver](scalar_t timestep, VariablesBasePtr vars) { io::HDF5File file("out_acc_carbo.hdf5", io::HDF5_OpenMode::OpenReadWrite); saver.save_timepoint(file, std::to_string(timestep), "/"); }; //solver::ReactiveTransportReturnCode retcode = react_solver.solve_timestep(0.01, vars); //retcode = react_solver.solve_timestep(1.0, vars); //std::cout << (int) retcode << std::endl; // retcode = react_solver.solve_timestep(5.0, vars); // std::cout << (int) retcode << std::endl; // retcode = react_solver.solve_timestep(10.0, vars); // std::cout << (int) retcode << std::endl; // retcode = react_solver.solve_timestep(10.0, vars); // std::cout << (int) retcode << std::endl; solver::ReactiveTransportRunner runner(react_solver, 0.001, 10, simul_info); runner.set_output_policy(out_pol); runner.run_until(30*24*3600, vars); std::cout << vars->get_liquid_saturation().variable << std::endl; std::cout << vars->get_liquid_saturation()(1) - vars->get_liquid_saturation()(2) << std::endl; std::cout << vars->get_pressure_main_variables(id_co2).variable << std::endl; } void UpscalingStagger::initialize(VariablesBasePtr var) { UnsaturatedVariables* true_var = static_cast(var.get()); // set relative models true_var->set_capillary_pressure_model( [this](index_t x, scalar_t sat){ return this->capillary_pressure(x, sat); }); true_var->set_relative_liquid_permeability_model( [this](index_t x, scalar_t sat){ return this->relative_liquid_permeability(x, sat); }); true_var->set_relative_liquid_diffusivity_model( [this](index_t x, scalar_t sat){ return this->relative_liquid_diffusivity(x, sat); }); true_var->set_relative_gas_diffusivity_model( [this](index_t x, scalar_t sat){ return this->relative_gas_diffusivity(x, sat); }); true_var->set_vapor_pressure_model( [this](index_t x, scalar_t sat){ return this->vapor_pressure(x, sat); }); // initialization SecondaryTransientVariable& porosity = true_var->get_porosity(); SecondaryVariable& permeability = true_var->get_liquid_permeability(); SecondaryVariable& liq_diffusivity = true_var->get_liquid_diffusivity(); SecondaryVariable& gas_diffusivity = true_var->get_resistance_gas_diffusivity(); true_var->set_relative_variables(); gas_diffusivity(0) = resistance_gas_diffusivity(1.0); true_var->get_relative_gas_diffusivity()(0) = relative_gas_diffusivity(0, 1.0); for (index_t node=1; nodeget_mesh()->nb_nodes(); ++node) { scalar_t phi = porosity(node); permeability(node) = intrinsic_permeability(phi); liq_diffusivity(node) = intrinsic_liquid_diffusivity(phi); gas_diffusivity(node) = resistance_gas_diffusivity(phi); //true_var->set_relative_variables(node); } } scalar_t UpscalingStagger::capillary_pressure(index_t node, scalar_t saturation) { if (saturation == 0.0) return 0.0; - scalar_t coeff = - constants::water_density_25 - * constants::gas_constant - * units::celsius(25.0) - / 18.015e-3; - return coeff*ln_hr(saturation); + return laws::capillary_pressure_BaroghelBouny(saturation, m_a, m_b); } scalar_t UpscalingStagger::vapor_pressure(index_t node, scalar_t saturation) { scalar_t pv; if (saturation == 0.0) pv = 0.0; else if (saturation >= 1.0) pv = m_pv_sat; - else pv = m_pv_sat*std::exp(ln_hr(saturation)); + else pv = m_pv_sat*laws::invert_kelvin_equation( + capillary_pressure(node, saturation)); return pv; } scalar_t UpscalingStagger::relative_liquid_diffusivity( index_t node, scalar_t saturation ) { if (saturation == 0.0) return 0.0; return std::pow(saturation, m_exp_r_l_d); } scalar_t UpscalingStagger::relative_liquid_permeability( index_t node, scalar_t saturation ) { if (saturation == 0.0) return 0.0; - scalar_t tmp = std::pow(saturation, 1.0/m_exp_r_l_k); - scalar_t tmp2 = std::pow(1.0 - tmp, m_exp_r_l_k); - tmp = std::pow(1.0 - tmp2, 2); - tmp2 = tmp*std::sqrt(saturation); - return tmp2; + return laws::relative_liquid_permeability_Mualem(saturation, 1/m_b); } scalar_t UpscalingStagger::relative_gas_diffusivity(index_t node, scalar_t saturation) { if (saturation == 0.0) return 1.0; return std::pow(1.0 - saturation, m_exp_r_g_d); } -scalar_t UpscalingStagger::ln_hr(scalar_t saturation) -{ - scalar_t tmp = std::pow(saturation, -1.0/m_beta) - 1.0; - scalar_t tmp2 = std::pow(tmp, m_beta/(1.0-m_beta)); - tmp = -1.0/m_alpha*tmp2; - return tmp; -} - scalar_t UpscalingStagger::intrinsic_permeability(scalar_t porosity) { scalar_t tmp_1 = porosity / m_phi_0; scalar_t tmp_2 = (1.0 - m_phi_0) / (1.0 - porosity); tmp_1 = std::pow(tmp_1, 3); tmp_2 = std::pow(tmp_2, 2); return m_k_0 * tmp_1 * tmp_2; } scalar_t UpscalingStagger::intrinsic_liquid_diffusivity(scalar_t porosity) { return m_d_0 * std::exp(9.95*porosity); } scalar_t UpscalingStagger::resistance_gas_diffusivity(scalar_t porosity) { return std::pow(porosity, m_exp_i_g_d); } void UpscalingStagger::initialize_timestep(scalar_t dt, VariablesBasePtr var) { } StaggerReturnCode UpscalingStagger::restart_timestep(VariablesBasePtr var) { UnsaturatedVariables* true_var = static_cast(var.get()); SecondaryTransientVariable& porosity = true_var->get_porosity(); SecondaryVariable& permeability = true_var->get_liquid_permeability(); SecondaryVariable& liq_diffusivity = true_var->get_liquid_diffusivity(); SecondaryVariable& gas_diffusivity = true_var->get_resistance_gas_diffusivity(); for (index_t node=1; nodeget_mesh()->nb_nodes(); ++node) { scalar_t phi = porosity(node); permeability(node) = intrinsic_permeability(phi); liq_diffusivity(node) = intrinsic_liquid_diffusivity(phi); gas_diffusivity(node) = resistance_gas_diffusivity(phi); true_var->set_relative_variables(node); } return StaggerReturnCode::ResidualMinimized; } specmicp::RawDatabasePtr get_database() { specmicp::database::Database thedatabase(CEMDATA_PATH); std::map swap = {{"HCO3[-]", "CO2"},}; thedatabase.swap_components(swap); database::RawDatabasePtr raw_data = thedatabase.get_database(); return raw_data; } AdimensionalSystemSolution get_initial_condition( database::RawDatabasePtr the_database, const units::UnitsSet& units_set ) { Formulation initial_formulation; scalar_t scaling = 0.9e-3; - initial_formulation.set_mass_solution(0.6*scaling); + initial_formulation.set_mass_solution(0.5*scaling); initial_formulation.add_mineral("C3S", 5.5*scaling); initial_formulation.add_mineral("C2S", 1.25*scaling); initial_formulation.add_mineral("Calcite", 0.05*scaling); Dissolver dissolver(the_database); UpscalingStagger upsstag_init; // just for vapor pressure model AdimensionalSystemConstraints constraints; constraints.set_total_concentrations( dissolver.dissolve(initial_formulation, true) ); AdimensionalSystemSolver the_solver(the_database, constraints); the_solver.get_options().units_set = units_set; the_solver.get_options().system_options.non_ideality_tolerance = 1e-12; micpsolver::MiCPSolverOptions* solver_options = &the_solver.get_options().solver_options; //solver_options.maxstep = 10.0; solver_options->maxiter_maxstep = 200; solver_options->set_tolerance(1e-9); Vector x; the_solver.initialise_variables(x, 0.8, -4); micpsolver::MiCPPerformance perf = the_solver.solve(x); if (perf.return_code < micpsolver::MiCPSolverReturnCode::Success) { throw std::runtime_error("Error : failed to solve initial conditions"); } constraints.set_water_partial_pressure_model( std::bind(&UpscalingStagger::vapor_pressure, &upsstag_init, 1, std::placeholders::_1) ); the_solver = AdimensionalSystemSolver (the_database, constraints); the_solver.get_options().units_set = units_set; the_solver.get_options().system_options.non_ideality_tolerance = 1e-12; solver_options = &the_solver.get_options().solver_options; //solver_options.maxstep = 10.0; solver_options->maxiter_maxstep = 200; solver_options->set_tolerance(1e-9); perf = the_solver.solve(x); if (perf.return_code < micpsolver::MiCPSolverReturnCode::Success) { throw std::runtime_error("Error : failed to solve initial conditions"); } return the_solver.get_raw_solution(x); } diff --git a/src/physics/constants.hpp b/src/physics/constants.hpp index ddab2fb..802470f 100644 --- a/src/physics/constants.hpp +++ b/src/physics/constants.hpp @@ -1,59 +1,61 @@ /*------------------------------------------------------------------------------- Copyright (c) 2014,2015 F. Georget , Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -----------------------------------------------------------------------------*/ #ifndef SPECMICP_CONSTANTS_CONSTANTS_HPP #define SPECMICP_CONSTANTS_CONSTANTS_HPP //! \file constants.hpp //! \brief Contains physical constants namespace specmicp { //! \namespace specmicp::constants //! \brief Constants used by SpecMiCP namespace constants { constexpr double gas_constant = 8.3144621; //!< Perfect gas constant (J / K mol) constexpr double water_density_25 = 997.05; //!< Density of water (kg / m^3) constexpr double water_density_25_kgl = 0.99705; //!< Density of water (kg / L) constexpr double water_viscosity = 8.937e-4; //!< Viscosity water (Pa.s) +constexpr double water_molar_mass = 18.01528e-3; //!< Molar mass water (kg / mol) + constexpr double Adebye = 0.5092; //!< constant 'A' for Debye-Huckel activity coefficient constexpr double Bdebye = 0.3283; //!< constant 'B' for Debye-Huckel activity coefficient constexpr double faraday_constant = 9.6485339924e4 ; //! the Faraday constant (C / mol) } // end namespace constants } // end namespace specmicp #endif // SPECMICP_CONSTANTS_CONSTANTS_HPP diff --git a/src/physics/constants.hpp b/src/physics/unsaturated_laws.cpp similarity index 53% copy from src/physics/constants.hpp copy to src/physics/unsaturated_laws.cpp index ddab2fb..7857c21 100644 --- a/src/physics/constants.hpp +++ b/src/physics/unsaturated_laws.cpp @@ -1,59 +1,89 @@ -/*------------------------------------------------------------------------------- +/*------------------------------------------------------------------------------ -Copyright (c) 2014,2015 F. Georget , Princeton University +Copyright (c)2015 F. Georget , Princeton University All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ------------------------------------------------------------------------------*/ +------------------------------------------------------------------------------*/ -#ifndef SPECMICP_CONSTANTS_CONSTANTS_HPP -#define SPECMICP_CONSTANTS_CONSTANTS_HPP +#include "unsaturated_laws.hpp" -//! \file constants.hpp -//! \brief Contains physical constants +#include "laws.hpp" +#include "units.hpp" namespace specmicp { -//! \namespace specmicp::constants -//! \brief Constants used by SpecMiCP -namespace constants { - -constexpr double gas_constant = 8.3144621; //!< Perfect gas constant (J / K mol) - -constexpr double water_density_25 = 997.05; //!< Density of water (kg / m^3) -constexpr double water_density_25_kgl = 0.99705; //!< Density of water (kg / L) - -constexpr double water_viscosity = 8.937e-4; //!< Viscosity water (Pa.s) - -constexpr double Adebye = 0.5092; //!< constant 'A' for Debye-Huckel activity coefficient -constexpr double Bdebye = 0.3283; //!< constant 'B' for Debye-Huckel activity coefficient - -constexpr double faraday_constant = 9.6485339924e4 ; //! the Faraday constant (C / mol) - -} // end namespace constants -} // end namespace specmicp - -#endif // SPECMICP_CONSTANTS_CONSTANTS_HPP +namespace laws { + +scalar_t kelvin_equation(scalar_t rh) +{ + return -(constants::gas_constant * units::celsius(25.0) * + constants::water_density_25 / constants::water_molar_mass + )*std::log(rh); +} + +scalar_t invert_kelvin_equation(scalar_t pc) +{ + return std::exp( -constants::water_molar_mass / + ( constants::gas_constant * + units::celsius(25.0) * + constants::water_density_25 + ) * pc + ); +} + +scalar_t capillary_pressure_BaroghelBouny( + scalar_t saturation, + scalar_t a, + scalar_t b + ) +{ + auto tmp = std::pow(saturation, -b) - 1.0; + return a*std::pow(tmp, 1.0-1.0/b); +} + + +scalar_t relative_liquid_permeability_Mualem( + scalar_t saturation, + scalar_t m + ) +{ + auto tmp = 1.0 - std::pow(saturation, 1.0/m); + auto tmp2 = 1.0 - std::pow(tmp, m); + return std::sqrt(saturation) * std::pow(tmp2, 2); +} + +scalar_t relative_gas_permeability_Mualem( + scalar_t saturation, + scalar_t m + ) +{ + auto tmp = 1.0 - std::pow(saturation, 1.0/m); + return std::sqrt(1.0-saturation)*std::pow(tmp, 2*m); +} + +} //end namespace laws +} //end namespace specmicp diff --git a/src/physics/unsaturated_laws.hpp b/src/physics/unsaturated_laws.hpp new file mode 100644 index 0000000..8020866 --- /dev/null +++ b/src/physics/unsaturated_laws.hpp @@ -0,0 +1,115 @@ +/*------------------------------------------------------------------------------ + +Copyright (c)2015 F. Georget , Princeton University +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +are permitted provided that the following conditions are met: + +1. Redistributions of source code must retain the above copyright notice, this +list of conditions and the following disclaimer. + +2. Redistributions in binary form must reproduce the above copyright notice, +this list of conditions and the following disclaimer in the documentation and/or +other materials provided with the distribution. + +3. Neither the name of the copyright holder nor the names of its contributors +may be used to endorse or promote products derived from this software without +specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR +ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON +ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +------------------------------------------------------------------------------*/ + +#ifndef SPECMICP_PHYSICS_UNSATURATEDLAWS_HPP +#define SPECMICP_PHYSICS_UNSATURATEDLAWS_HPP + +//! \file physics/unsaturated_laws.hpp +//! \brief Common laws and models in unsaturated porous medium + +#include "../types.hpp" + +namespace specmicp { +namespace laws { + +//! \brief Return the capillary pressure as function of the relative humidity +scalar_t SPECMICP_DLL_PUBLIC kelvin_equation(scalar_t rh); + +//! \brief Return the relative humidity as function of the capillary pressure +scalar_t SPECMICP_DLL_PUBLIC invert_kelvin_equation(scalar_t pc); + +//! \brief Capillary Pressure model used by V. Baroghel-Bouny et al. +//! +//! \f$ Pc(S_l) = a (S_l^{-b} - 1)^{1-1/b} \f$ +//! +//! \param saturation the liquid saturation +//! \param a fitting coefficient, unit of pressure +//! \param b fitting coefficient, dimensioneless, >1 +scalar_t SPECMICP_DLL_PUBLIC capillary_pressure_BaroghelBouny( + scalar_t saturation, + scalar_t a, + scalar_t b + ); + +//! \brief The Mualem model for relative liquid permeability +//! +//! \f$ k_{rl}(S_l) = \sqrt{S_l}\left(1 - \left(1-S_l^{1/m}\right)^m\right)^2 \f$ +//! +//! \param saturation the liquid saturation +//! \param m a fitting coefficient +//! +//! In Baroghel-Bouny et al. m=1.0/b +scalar_t SPECMICP_DLL_PUBLIC relative_liquid_permeability_Mualem( + scalar_t saturation, + scalar_t m + ); + +//! \brief The mualem Model for relative gas permeability +//! +//! \f$ k_{rg} = \sqrt{1-S_l} \left(1 - S)l^{1/m} \right)^{2m} +scalar_t SPECMICP_DLL_PUBLIC relative_gas_permeability_Mualem( + scalar_t saturation, + scalar_t m + ); + +//! \brief Type of a model of the saturation +using saturation_model_f = std::function; + +//! \brief Return a model function of the saturation +//! +//! \tparam F a model where saturation is the first parameter +//! \tparam ...Args Parameter pack representing the parameters of the model +//! \return a function taking only the saturation +//! +//! Bind the argument of the model +template +saturation_model_f make_saturation_model(F func, Args...args) +{ + return [func,args...](scalar_t saturation) -> scalar_t { + return func(saturation, args...); + }; +} + +//! \brief Return a relative humidity model from a capillary pressure model +//! +//! Use the kelvin equation +inline saturation_model_f make_rh_model(saturation_model_f cap_pressure) +{ + return [cap_pressure](scalar_t saturation) -> scalar_t { + return invert_kelvin_equation(cap_pressure(saturation)); + }; +} + +} //end namespace laws +} //end namespace specmicp + +#endif // SPECMICP_PHYSICS_UNSATURATEDLAWS_HPP diff --git a/src/utils/CMakeLists.txt b/src/utils/CMakeLists.txt index a49ab41..7e0ec8b 100644 --- a/src/utils/CMakeLists.txt +++ b/src/utils/CMakeLists.txt @@ -1,164 +1,165 @@ # headers only module # ------------------- add_custom_target(utils_inc SOURCES #log.hpp # sparse solvers # -------------- sparse_solvers/sparse_solver_base.hpp sparse_solvers/sparse_solver.hpp sparse_solvers/sparse_solver_structs.hpp sparse_solvers/sparse_qr.hpp sparse_solvers/sparse_lu.hpp sparse_solvers/sparse_bicgstab.hpp sparse_solvers/sparse_gmres.hpp range_iterator.hpp range_iterator.inl options_handler.hpp perfs_handler.hpp compat.hpp value_checker.hpp scope_guard.hpp plugins/plugin_api_version.h plugins/plugin_interface.h plugins/plugin_types.hpp ) set(SPECMICP_COMMON_LIBRARY_FILES log.cpp dateandtime.cpp timer.cpp moving_average.cpp plugins/dynamic_library.cpp plugins/plugin_base.cpp plugins/module_base.cpp plugins/plugin_manager.cpp ../physics/laws.cpp + ../physics/unsaturated_laws.cpp ../physics/units.cpp ../physics/io/units.cpp io/format.cpp io/csv_formatter.cpp io/yaml.cpp cli/parser.cpp ) # hdf5 # ----- if (HDF5_FOUND) list(APPEND SPECMICP_COMMON_LIBRARY_FILES io/specmicp_hdf5.cpp ) include_directories(${HDF5_INCLUDE_DIRS}) set_source_files_properties(io/specmicp_hdf5.cpp PROPERTIES COMPILE_DEFINITIONS HDF5_DEFINITIONS) add_custom_target(utils_io_hdf5_inc SOURCES io/hdf5_eigen.hpp io/hdf5_eigen.inl ) endif() set_pgo_flag(${SPECMICP_COMMON_LIBRARY_FILES}) add_library(objspecmicp_common OBJECT ${SPECMICP_COMMON_LIBRARY_FILES}) set_property(TARGET objspecmicp_common PROPERTY POSITION_INDEPENDENT_CODE 1) add_library(specmicp_common SHARED $) if (UNIX) set(SPECMICP_COMMON_LINK dl) else() message(FATAL_ERROR "Plugin system only for POSIX at this time !") endif() if (HDF5_FOUND) set(SPECMICP_COMMON_LINK "${HDF5_LIBRARIES};${SPECMICP_COMMON_LINK}") endif() target_link_libraries(specmicp_common ${SPECMICP_COMMON_LINK}) install(TARGETS specmicp_common LIBRARY DESTINATION ${LIBRARY_INSTALL_DIR} ) set(UTILS_INCLUDE_LIST range_iterator.hpp range_iterator.inl log.hpp options_handler.hpp perfs_handler.hpp moving_average.hpp timer.hpp dateandtime.hpp compat.hpp scope_guard.hpp ) set(UTILS_SPARSE_INCLUDE_LIST # sparse solvers # -------------- sparse_solvers/sparse_solver.hpp sparse_solvers/sparse_solver_base.hpp sparse_solvers/sparse_solver_structs.hpp sparse_solvers/sparse_qr.hpp sparse_solvers/sparse_lu.hpp sparse_solvers/sparse_bicgstab.hpp sparse_solvers/sparse_gmres.hpp ) set(UTILS_IO_INCLUDE_LIST io/format.hpp io/csv_formatter.hpp io/yaml.hpp ) set(UTILS_CLI_INCLUDE_LIST cli/parser.hpp ) install(FILES ${UTILS_INCLUDE_LIST} DESTINATION ${INCLUDE_INSTALL_DIR}/utils ) install(FILES ${UTILS_SPARSE_INCLUDE_LIST} DESTINATION ${INCLUDE_INSTALL_DIR}/utils/sparse_solvers ) install(FILES ${UTILS_IO_INCLUDE_LIST} DESTINATION ${INCLUDE_INSTALL_DIR}/utils/io ) # static libraries # ---------------- if(SPECMICP_BUILD_STATIC) add_library(specmicp_common_static STATIC $ ) install(TARGETS specmicp_common_static ARCHIVE DESTINATION ${STATIC_LIBRARY_INSTALL_DIR} ) else() add_library(specmicp_common_static EXCLUDE_FROM_ALL STATIC $ ) endif() target_link_libraries(specmicp_common_static ${SPECMICP_COMMON_LINK}) set_target_properties(specmicp_common_static PROPERTIES OUTPUT_NAME specmicp_common) diff --git a/tests/common/laws.cpp b/tests/common/laws.cpp index e355767..5c6fa53 100644 --- a/tests/common/laws.cpp +++ b/tests/common/laws.cpp @@ -1,56 +1,114 @@ #include "catch.hpp" #include "physics/laws.hpp" +#include "physics/unsaturated_laws.hpp" + +#include using namespace specmicp; TEST_CASE("Gas perfect") { SECTION("Pressure") { REQUIRE(laws::pressure_perfect_gas(1.0, 1.0, 1.0) == Approx(constants::gas_constant)); REQUIRE(laws::pressure_perfect_gas(2.0, 1.0, 1.0) == Approx(2.0*constants::gas_constant)); REQUIRE(laws::pressure_perfect_gas(1.0, 3.0, 1.0) == Approx(constants::gas_constant/3.0)); REQUIRE(laws::pressure_perfect_gas(1.0, 1.0, 4.0) == Approx(constants::gas_constant*4.0)); REQUIRE(laws::pressure_perfect_gas(2.0, 3.0, 4.0) == Approx(2.0*constants::gas_constant*4.0/3.0)); } SECTION("Mole") { REQUIRE(laws::mole_perfect_gas(1.0, 1.0, 1.0) == Approx(1.0/constants::gas_constant)); REQUIRE(laws::mole_perfect_gas(2.0, 1.0, 1.0) == Approx(2.0/constants::gas_constant)); REQUIRE(laws::mole_perfect_gas(1.0, 3.0, 1.0) == Approx(3.0/constants::gas_constant)); REQUIRE(laws::mole_perfect_gas(1.0, 1.0, 4.0) == Approx(1.0/(4.0*constants::gas_constant))); REQUIRE(laws::mole_perfect_gas(2.0, 3.0, 4.0) == Approx(2.0*3.0/(4.0*constants::gas_constant))); } } TEST_CASE("Debye-Huckel") { SECTION("Debye-Huckel") { REQUIRE(laws::debye_huckel(1.0, 1.0, 1.0) == Approx(-(constants::Adebye)/(1+constants::Bdebye))); REQUIRE(laws::debye_huckel(1.0, 2.0, 1.0) == Approx(-(constants::Adebye*4)/(1+constants::Bdebye))); REQUIRE(laws::debye_huckel(0.1, 1.0, 1.0) == Approx(-(constants::Adebye*0.1)/(1+constants::Bdebye*0.1))); REQUIRE(laws::debye_huckel(0.1, 2.0, 1.0) == Approx(-(constants::Adebye*0.1*4)/(1+constants::Bdebye*0.1))); REQUIRE(laws::debye_huckel(1.0, 1.0, 2.0) == Approx(-(constants::Adebye)/(1+constants::Bdebye*1.0*2.0))); REQUIRE(laws::debye_huckel(0.1, 2.0, 3.0) == Approx(-(constants::Adebye*0.1*4)/(1+constants::Bdebye*0.1*3.0))); } SECTION("Extended Debye-Huckel") { REQUIRE(laws::extended_debye_huckel(1.0, 1.0, 1.0, 1.0, 1.0) == Approx(-(constants::Adebye)/(1.0+constants::Bdebye)+1.0*1.0)); REQUIRE(laws::extended_debye_huckel(1.0, 1.0, 2.0, 1.0, 1.0) == Approx(-(constants::Adebye*4.0)/(1.0+constants::Bdebye)+1.0*1.0)); REQUIRE(laws::extended_debye_huckel(1.0, 1.0, 1.0, 1.0, 2.0) == Approx(-(constants::Adebye*1.0)/(1.0+constants::Bdebye)+1.0*2.0)); REQUIRE(laws::extended_debye_huckel(0.1, 0.01, 1.0, 1.0, 2.0) == Approx(-(constants::Adebye*1.0*0.01)/(1.0+constants::Bdebye*0.01)+2.0*0.1)); REQUIRE(laws::extended_debye_huckel(0.01, 0.1, -2.0, 3.0, 2.0) == Approx(-(4*constants::Adebye*1.0*0.1)/(1.0+constants::Bdebye*0.1*3.0)+2.0*0.01)); REQUIRE(laws::extended_debye_huckel(0.01, -2.0, 3.0, 2.0) == Approx(-(4*constants::Adebye*1.0*0.1)/(1.0+constants::Bdebye*0.1*3.0)+2.0*0.01)); } } + +TEST_CASE("Unsaturated laws", "[laws],[unsaturated],[physics]") { + SECTION("make_saturation_model") { + auto identity = [](scalar_t saturation) -> scalar_t {return saturation;}; + auto model = laws::make_saturation_model(identity); + CHECK(model(1.0) == 1.0); + CHECK(model(0.5) == 0.5); + + auto mult = [](scalar_t sat, scalar_t mult) -> scalar_t { + return mult*sat; + }; + auto model2 = laws::make_saturation_model(mult, 2.0); + CHECK(model2(1.0) == Approx(2.0)); + CHECK(model2(0.5) == Approx(1.0)); + } + + SECTION("Baroghel Bouny model") { + auto cap = laws::make_saturation_model( + &laws::capillary_pressure_BaroghelBouny, + 1.0e8, 2.0); + CHECK(cap(1.0) == 0.0); + CHECK(cap(0.5) == Approx(1.732050808e8)); + + auto rh = laws::make_rh_model(cap); + + CHECK(rh(1.0) == 1.0); + CHECK(rh(0.0) == 0.0); + CHECK(rh(0.5) == Approx(0.2832770066).epsilon(1e-3)); + + cap = laws::make_saturation_model( + &laws::capillary_pressure_BaroghelBouny, + 2.0, 2.0); + CHECK(cap(0.5) == Approx(2.0*1.732050808)); + + } + + SECTION("Mualem model") { + auto krl = laws::make_saturation_model( + &laws::relative_liquid_permeability_Mualem, + 2.0 + ); + CHECK(krl(1.0) == Approx(1.0)); + CHECK(krl(0.0) == Approx(0.0)); + CHECK(krl(0.5) == Approx(0.5909902577)); + + auto krg = laws::make_saturation_model( + &laws::relative_gas_permeability_Mualem, + 2.0 + ); + + CHECK(krg(1.0) == Approx(0.0)); + CHECK(krg(0.0) == Approx(1.0)); + CHECK(krg(0.5) == Approx(0.005203820043)); + } +}