diff --git a/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp b/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp index c688a1211563ed077621eae0afc579d3dc4d08d5..8652773a330145eb3b47c8952624f78c7f304589 100644 --- a/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp +++ b/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp @@ -112,22 +112,6 @@ std::unique_ptr<Process> createHydroMechanicsProcess( MaterialLib::Solids::createConstitutiveRelations<DisplacementDim>( parameters, local_coordinate_system, config); - // Intrinsic permeability - auto& intrinsic_permeability = ParameterLib::findParameter<double>( - config, - //! \ogs_file_param_special{prj__processes__process__HYDRO_MECHANICS__intrinsic_permeability} - "intrinsic_permeability", parameters, 1, &mesh); - - DBUG("Use '%s' as intrinsic conductivity parameter.", - intrinsic_permeability.name.c_str()); - - // Fluid density - auto& fluid_density = ParameterLib::findParameter<double>( - config, - //! \ogs_file_param_special{prj__processes__process__HYDRO_MECHANICS__fluid_density} - "fluid_density", parameters, 1, &mesh); - DBUG("Use '%s' as fluid density parameter.", fluid_density.name.c_str()); - // Specific body force Eigen::Matrix<double, DisplacementDim, 1> specific_body_force; { @@ -150,10 +134,12 @@ std::unique_ptr<Process> createHydroMechanicsProcess( auto media_map = MaterialPropertyLib::createMaterialSpatialDistributionMap(media, mesh); - std::array const requiredGasProperties = {MaterialPropertyLib::viscosity}; + std::array const requiredGasProperties = { + MaterialPropertyLib::viscosity, MaterialPropertyLib::density}; std::array const requiredSolidProperties = { MaterialPropertyLib::porosity, MaterialPropertyLib::biot_coefficient, - MaterialPropertyLib::density}; + MaterialPropertyLib::density, MaterialPropertyLib::permeability}; + for (auto const& m : media) { checkRequiredProperties(m.second->phase("Gas"), requiredGasProperties); @@ -206,8 +192,7 @@ std::unique_ptr<Process> createHydroMechanicsProcess( HydroMechanicsProcessData<DisplacementDim> process_data{ materialIDs(mesh), std::move(media_map), std::move(solid_constitutive_relations), - initial_stress, intrinsic_permeability, - fluid_density, specific_body_force, + initial_stress, specific_body_force, fluid_compressibility, reference_temperature, specific_gas_constant, fluid_type}; diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h index 7bc079a069a7115ffb234c0e8f1a647418c7c75e..58d3263baf347ba99e1691ffba4ca7bda57756df 100644 --- a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h +++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h @@ -215,11 +215,12 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement, auto& eps = _ip_data[ip].eps; auto const& sigma_eff = _ip_data[ip].sigma_eff; + auto const K = solid_phase.property(MPL::PropertyType::permeability) + .template value<double>(vars, x_position, t); auto const mu = gas_phase.property(MPL::PropertyType::viscosity) .template value<double>(vars, x_position, t); - double const K_over_mu = - _process_data.intrinsic_permeability(t, x_position)[0] / mu; + auto const K_over_mu = K / mu; auto const alpha = solid_phase.property(MPL::PropertyType::biot_coefficient) @@ -232,9 +233,9 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement, (_process_data.fluid_type == FluidType::Fluid_Type::IDEAL_GAS) ? N_p.dot(p) : std::numeric_limits<double>::quiet_NaN(); - double const rho_fr = - _process_data.getFluidDensity(t, x_position, p_fr); double const beta_p = _process_data.getFluidCompressibility(p_fr); + auto const rho_fr = gas_phase.property(MPL::PropertyType::density) + .template value<double>(vars, x_position, t); auto const porosity = solid_phase.property(MPL::PropertyType::porosity) .template value<double>(vars, x_position, t); auto const& identity2 = MathLib::KelvinVector::Invariants< @@ -347,19 +348,24 @@ HydroMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure, auto const& medium = _process_data.media_map->getMedium(_element.getID()); auto const& gas_phase = medium->phase("Gas"); + auto const& solid_phase = medium->phase("Solid"); MPL::VariableArray vars; for (unsigned ip = 0; ip < n_integration_points; ip++) { x_position.setIntegrationPoint(ip); + auto const K = solid_phase.property(MPL::PropertyType::permeability) + .template value<double>(vars, x_position, t); + auto const mu = gas_phase.property(MPL::PropertyType::viscosity) .template value<double>(vars, x_position, t); - double const K_over_mu = - _process_data.intrinsic_permeability(t, x_position)[0] / mu; + auto const K_over_mu = K / mu; + + auto const rho_fr = gas.property(MPL::PropertyType::density) + .template value<double>(vars, x_position, t); - auto const rho_fr = _process_data.fluid_density(t, x_position)[0]; auto const& b = _process_data.specific_body_force; // Compute the velocity @@ -434,11 +440,13 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement, auto const& N_p = _ip_data[ip].N_p; auto const& dNdx_p = _ip_data[ip].dNdx_p; + auto const K = solid_phase.property(MPL::PropertyType::permeability) + .template value<double>(vars, x_position, t); auto const mu = gas_phase.property(MPL::PropertyType::viscosity) .template value<double>(vars, x_position, t); - double const K_over_mu = - _process_data.intrinsic_permeability(t, x_position)[0] / mu; + auto const K_over_mu = K / mu; + auto const alpha_b = solid_phase.property(MPL::PropertyType::biot_coefficient) .template value<double>(vars, x_position, t); @@ -447,9 +455,9 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement, (_process_data.fluid_type == FluidType::Fluid_Type::IDEAL_GAS) ? N_p.dot(p) : std::numeric_limits<double>::quiet_NaN(); - double const rho_fr = - _process_data.getFluidDensity(t, x_position, p_fr); double const beta_p = _process_data.getFluidCompressibility(p_fr); + auto const rho_fr = gas_phase.property(MPL::PropertyType::density) + .template value<double>(vars, x_position, t); auto const porosity = solid_phase.property(MPL::PropertyType::porosity) .template value<double>(vars, x_position, t); auto const K_S = solid_material.getBulkModulus(t, x_position); @@ -515,6 +523,7 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement, x_position.setElementID(_element.getID()); auto const& medium = _process_data.media_map->getMedium(_element.getID()); + auto const& gas_phase = medium->phase("Gas"); auto const& solid_phase = medium->phase("Solid"); MPL::VariableArray vars; @@ -547,9 +556,11 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement, auto const alpha = solid_phase.property(MPL::PropertyType::biot_coefficient) .template value<double>(vars, x_position, t); + auto const rho_sr = solid_phase.property(MPL::PropertyType::density) .template value<double>(vars, x_position, t); - auto const rho_fr = _process_data.fluid_density(t, x_position)[0]; + auto const rho_fr = gas_phase.property(MPL::PropertyType::density) + .template value<double>(vars, x_position, t); auto const porosity = solid_phase.property(MPL::PropertyType::porosity) .template value<double>(vars, x_position, t); auto const& b = _process_data.specific_body_force; diff --git a/ProcessLib/HydroMechanics/HydroMechanicsProcessData.h b/ProcessLib/HydroMechanics/HydroMechanicsProcessData.h index 74b694ec7c15eafdbd8a5b01b13a34ee8acca9d3..f9cd065565a8db1ee13428c61e6310c72c973d74 100644 --- a/ProcessLib/HydroMechanics/HydroMechanicsProcessData.h +++ b/ProcessLib/HydroMechanics/HydroMechanicsProcessData.h @@ -51,11 +51,6 @@ struct HydroMechanicsProcessData /// representation of length 4 or 6, ParameterLib::Parameter<double>. ParameterLib::Parameter<double> const* const initial_stress; - /// Permeability of the solid. A scalar quantity, - /// ParameterLib::Parameter<double>. - ParameterLib::Parameter<double> const& intrinsic_permeability; - /// Fluid's density. A scalar quantity, ParameterLib::Parameter<double>. - ParameterLib::Parameter<double> const& fluid_density; /// Specific body forces applied to solid and fluid. /// It is usually used to apply gravitational forces. /// A vector of displacement dimension's length. @@ -80,22 +75,6 @@ struct HydroMechanicsProcessData /// incompressible_fluid, compressible_fluid, ideal_gas FluidType::Fluid_Type const fluid_type; - /// will be removed after linking with MPL - double getFluidDensity( - double const& t, ParameterLib::SpatialPosition const& x_position, - double const& p_fr) - { - if (fluid_type == FluidType::Fluid_Type::INCOMPRESSIBLE_FLUID || - fluid_type == FluidType::Fluid_Type::COMPRESSIBLE_FLUID) - { - return fluid_density(t, x_position)[0]; - } - if (fluid_type == FluidType::Fluid_Type::IDEAL_GAS) - { - return p_fr / (specific_gas_constant * reference_temperature); - } - OGS_FATAL("unknown fluid type %d", static_cast<int> (fluid_type)); - } /// will be removed after linking with MPL double getFluidCompressibility(double const& p_fr)