Skip to content

Commit

Permalink
Merge pull request #2544 from zhangning737/master
Browse files Browse the repository at this point in the history
 Added nonequilibrium_stress to ThermoMechanics process
  • Loading branch information
endJunction committed Jul 5, 2019
2 parents 7dceab0 + dc81fde commit 29ce4b5
Show file tree
Hide file tree
Showing 20 changed files with 830 additions and 6 deletions.
Original file line number Diff line number Diff line change
@@ -0,0 +1,3 @@
Non-equilibrium initial state variables which are not explicitly in equilibrium
with any boundary conditions. They are not considered in the weak form but
affect the constitutive behaviour of the materials.
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
Non-equilibrium stress.
16 changes: 15 additions & 1 deletion ProcessLib/ThermoMechanics/CreateThermoMechanicsProcess.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -159,12 +159,26 @@ std::unique_ptr<Process> createThermoMechanicsProcess(

std::copy_n(b.data(), b.size(), specific_body_force.data());
}
// Non-equilibrium variables
ParameterLib::Parameter<double> const* nonequilibrium_stress = nullptr;
const auto& nonequilibrium_state_variables_config =
//! \ogs_file_param{prj__processes__process__SMALL_DEFORMATION__nonequilibrium_state_variables}
config.getConfigSubtreeOptional("nonequilibrium_state_variables");
if (nonequilibrium_state_variables_config)
{
nonequilibrium_stress = &ParameterLib::findParameter<double>(
*nonequilibrium_state_variables_config,
//! \ogs_file_param_special{prj__processes__process__SMALL_DEFORMATION__nonequilibrium_state_variables__stress}
"stress", parameters,
MathLib::KelvinVector::KelvinVectorDimensions<
DisplacementDim>::value);
}

ThermoMechanicsProcessData<DisplacementDim> process_data{
materialIDs(mesh), std::move(solid_constitutive_relations),
reference_solid_density, linear_thermal_expansion_coefficient,
specific_heat_capacity, thermal_conductivity,
specific_body_force};
specific_body_force, nonequilibrium_stress};

SecondaryVariableCollection secondary_variables;

Expand Down
3 changes: 3 additions & 0 deletions ProcessLib/ThermoMechanics/Tests.cmake
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,9 @@ if (NOT OGS_USE_MPI)
OgsTest(PROJECTFILE ThermoMechanics/CreepBGRa/Verification/m2_2Dloadbt/m2_2Dloadbt.prj RUNTIME 64)
OgsTest(PROJECTFILE ThermoMechanics/CreepBGRa/Verification/m2_3Dload/m2_3Dload.prj RUNTIME 24)
OgsTest(PROJECTFILE ThermoMechanics/CreepBGRa/Verification/m2_3Dloadbt/m2_3Dloadbt.prj RUNTIME 67)
OgsTest(PROJECTFILE ThermoMechanics/InitialStates/into_initial_state.prj)
OgsTest(PROJECTFILE ThermoMechanics/InitialStates/equilibrium_restart.prj)
OgsTest(PROJECTFILE ThermoMechanics/InitialStates/non_equilibrium_initial_state.prj)
endif()

AddTest(
Expand Down
38 changes: 35 additions & 3 deletions ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,12 @@ struct IntegrationPointData final
// mechanical strain
typename BMatricesType::KelvinVectorType eps_m, eps_m_prev;

// Non-equilibrium initial stress.
typename BMatricesType::KelvinVectorType sigma_neq =
BMatricesType::KelvinVectorType::Zero(
MathLib::KelvinVector::KelvinVectorDimensions<
DisplacementDim>::value);

MaterialLib::Solids::MechanicsBase<DisplacementDim> const& solid_material;
std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
DisplacementDim>::MaterialStateVariables>
Expand Down Expand Up @@ -164,6 +170,29 @@ class ThermoMechanicsLocalAssembler
ip_data.dNdx = shape_matrices[ip].dNdx;

_secondary_data.N[ip] = shape_matrices[ip].N;

// Initialize current time step values
if (_process_data.nonequilibrium_stress)
{
// Computation of non-equilibrium stress.
x_position.setCoordinates(MathLib::Point3d(
interpolateCoordinates<ShapeFunction, ShapeMatricesType>(
e, ip_data.N)));
std::vector<double> sigma_neq_data =
(*_process_data.nonequilibrium_stress)(
std::numeric_limits<
double>::quiet_NaN() /* time independent */,
x_position);
ip_data.sigma_neq =
Eigen::Map<typename BMatricesType::KelvinVectorType>(
sigma_neq_data.data(),
MathLib::KelvinVector::KelvinVectorDimensions<
DisplacementDim>::value,
1);
}
// Initialization from non-equilibrium sigma, which is zero by
// default, or is set to some value.
ip_data.sigma = ip_data.sigma_neq;
}
}

Expand Down Expand Up @@ -274,6 +303,7 @@ class ThermoMechanicsLocalAssembler

auto& sigma = _ip_data[ip].sigma;
auto const& sigma_prev = _ip_data[ip].sigma_prev;
auto const& sigma_neq = _ip_data[ip].sigma_neq;

auto& eps = _ip_data[ip].eps;
auto const& eps_prev = _ip_data[ip].eps_prev;
Expand Down Expand Up @@ -348,8 +378,9 @@ class ThermoMechanicsLocalAssembler
auto const& b = _process_data.specific_body_force;
local_rhs
.template block<displacement_size, 1>(displacement_index, 0)
.noalias() -=
(B.transpose() * sigma - N_u.transpose() * rho_s * b) * w;
.noalias() -= (B.transpose() * (sigma - sigma_neq) -
N_u.transpose() * rho_s * b) *
w;

//
// displacement equation, temperature part
Expand All @@ -359,7 +390,8 @@ class ThermoMechanicsLocalAssembler
if (_ip_data[ip].solid_material.getConstitutiveModel() ==
MaterialLib::Solids::ConstitutiveModel::CreepBGRa)
{
auto const s = Invariants::deviatoric_projection * sigma;
auto const s =
Invariants::deviatoric_projection * (sigma - sigma_neq);
double const norm_s = Invariants::FrobeniusNorm(s);
const double creep_coefficient =
_ip_data[ip]
Expand Down
7 changes: 5 additions & 2 deletions ProcessLib/ThermoMechanics/ThermoMechanicsProcessData.h
Original file line number Diff line number Diff line change
Expand Up @@ -42,15 +42,17 @@ struct ThermoMechanicsProcessData
linear_thermal_expansion_coefficient_,
ParameterLib::Parameter<double> const& specific_heat_capacity_,
ParameterLib::Parameter<double> const& thermal_conductivity_,
Eigen::Matrix<double, DisplacementDim, 1> const& specific_body_force_)
Eigen::Matrix<double, DisplacementDim, 1> const& specific_body_force_,
ParameterLib::Parameter<double> const* const nonequilibrium_stress_)
: material_ids(material_ids_),
solid_materials{std::move(solid_materials_)},
reference_solid_density(reference_solid_density_),
linear_thermal_expansion_coefficient(
linear_thermal_expansion_coefficient_),
specific_heat_capacity(specific_heat_capacity_),
thermal_conductivity(thermal_conductivity_),
specific_body_force(specific_body_force_)
specific_body_force(specific_body_force_),
nonequilibrium_stress(nonequilibrium_stress_)
{
}

Expand Down Expand Up @@ -79,6 +81,7 @@ struct ThermoMechanicsProcessData
ParameterLib::Parameter<double> const&
thermal_conductivity; // TODO To be changed as a matrix type variable.
Eigen::Matrix<double, DisplacementDim, 1> const specific_body_force;
ParameterLib::Parameter<double> const* const nonequilibrium_stress;
double dt = 0;
double t = 0;

Expand Down
241 changes: 241 additions & 0 deletions Tests/Data/ThermoMechanics/InitialStates/equilibrium_restart.prj
Original file line number Diff line number Diff line change
@@ -0,0 +1,241 @@
<?xml version="1.0" encoding="ISO-8859-1"?>
<OpenGeoSysProject>
<meshes>
<mesh>soil_column_q8_with_ic.vtu</mesh>
<mesh>soil_column_bottom.vtu</mesh>
<mesh>soil_column_top.vtu</mesh>
<mesh>soil_column_left.vtu</mesh>
<mesh>soil_column_right.vtu</mesh>
</meshes>
<processes>
<process>
<name>TM</name>
<type>THERMO_MECHANICS</type>
<integration_order>3</integration_order>
<constitutive_relation>
<type>LinearElasticIsotropic</type>
<youngs_modulus>E</youngs_modulus>
<poissons_ratio>nu</poissons_ratio>
</constitutive_relation>
<reference_solid_density>rho_sr</reference_solid_density>
<linear_thermal_expansion_coefficient>alpha</linear_thermal_expansion_coefficient>
<specific_heat_capacity>cs</specific_heat_capacity>
<thermal_conductivity>lambda</thermal_conductivity>
<specific_body_force>0 -9.81</specific_body_force>
<process_variables>
<displacement>displacement</displacement>
<temperature>temperature</temperature>
</process_variables>
<secondary_variables>
<secondary_variable type="static" internal_name="sigma" output_name="sigma"/>
<secondary_variable type="static" internal_name="epsilon" output_name="epsilon"/>
</secondary_variables>
</process>
</processes>
<time_loop>
<processes>
<process ref="TM">
<nonlinear_solver>basic_newton</nonlinear_solver>
<convergence_criterion>
<type>DeltaX</type>
<norm_type>NORM2</norm_type>
<abstol>1e-12</abstol>
</convergence_criterion>
<time_discretization>
<type>BackwardEuler</type>
</time_discretization>
<time_stepping>
<type>FixedTimeStepping</type>
<t_initial>1</t_initial>
<t_end>2</t_end>
<timesteps>
<pair>
<repeat>2</repeat>
<delta_t>1</delta_t>
</pair>
</timesteps>
</time_stepping>
</process>
</processes>
<output>
<type>VTK</type>
<prefix>equilibrium_restart</prefix>
<timesteps>
<pair>
<repeat>2</repeat>
<each_steps>1</each_steps>
</pair>
</timesteps>
<variables>
<variable>displacement</variable>
<variable>temperature</variable>
<variable>sigma</variable>
<variable>epsilon</variable>
</variables>
</output>
</time_loop>
<parameters>
<parameter>
<name>E</name>
<type>Constant</type>
<value>1e9</value>
</parameter>
<parameter>
<name>nu</name>
<type>Constant</type>
<value>.0</value>
</parameter>
<parameter>
<name>rho_sr</name>
<type>Constant</type>
<value>2200.0</value>
</parameter>
<parameter>
<name>alpha</name>
<type>Constant</type>
<value>6e-6</value>
</parameter>
<parameter>
<name>cs</name>
<type>Constant</type>
<value>1</value>
</parameter>
<parameter>
<name>lambda</name>
<type>Constant</type>
<value>1</value>
</parameter>
<parameter>
<name>temperature_ic</name>
<type>Constant</type>
<value>0</value>
</parameter>
<parameter>
<name>displacement0</name>
<type>Constant</type>
<values>0 0</values>
</parameter>
<parameter>
<name>Dirichlet_left</name>
<type>Constant</type>
<value>0.</value>
</parameter>
<parameter>
<name>Zero</name>
<type>Constant</type>
<value>0.</value>
</parameter>
<parameter>
<name>maximum_top_pressure</name>
<type>Constant</type>
<value>-1e6</value>
</parameter>
<parameter>
<name>top_load</name>
<type>CurveScaled</type>
<curve>top_load_ramp</curve>
<parameter>maximum_top_pressure</parameter>
</parameter>
<parameter>
<mesh>soil_column_right</mesh>
<name>side_load</name>
<type>Function</type>
<expression>-0.8*2200.0*9.81*(100.0-y)</expression>
</parameter>
</parameters>
<curves>
<curve>
<name>top_load_ramp</name>
<coords>0.0 1.0 2.0</coords>
<values>0.0 0.5 1.0</values>
</curve>
</curves>
<process_variables>
<process_variable>
<name>displacement</name>
<components>2</components>
<order>2</order>
<initial_condition>displacement0</initial_condition>
<boundary_conditions>
<!-- fixed boundaries -->
<boundary_condition>
<mesh>soil_column_bottom</mesh>
<type>Dirichlet</type>
<component>1</component>
<parameter>Zero</parameter>
</boundary_condition>
<boundary_condition>
<mesh>soil_column_left</mesh>
<type>Dirichlet</type>
<component>0</component>
<parameter>Zero</parameter>
</boundary_condition>
<boundary_condition>
<mesh>soil_column_top</mesh>
<type>Neumann</type>
<component>1</component>
<parameter>top_load</parameter>
</boundary_condition>
<boundary_condition>
<mesh>soil_column_right</mesh>
<type>Neumann</type>
<component>0</component>
<parameter>side_load</parameter>
</boundary_condition>
</boundary_conditions>
</process_variable>
<process_variable>
<name>temperature</name>
<components>1</components>
<order>1</order>
<initial_condition>temperature_ic</initial_condition>
<boundary_conditions>
<boundary_condition>
<mesh>soil_column_left</mesh>
<type>Dirichlet</type>
<parameter>temperature_ic</parameter>
</boundary_condition>
</boundary_conditions>
</process_variable>
</process_variables>
<nonlinear_solvers>
<nonlinear_solver>
<name>basic_newton</name>
<type>Newton</type>
<max_iter>10</max_iter>
<linear_solver>general_linear_solver</linear_solver>
</nonlinear_solver>
</nonlinear_solvers>
<linear_solvers>
<linear_solver>
<name>general_linear_solver</name>
<lis>-i cg -p jacobi -tol 1e-16 -maxiter 10000</lis>
<eigen>
<solver_type>SparseLU</solver_type>
<precon_type>DIAGONAL</precon_type>
<!--solver_type>CG</solver_type>
<precon_type>DIAGONAL</precon_type>
<max_iteration_step>10000</max_iteration_step>
<error_tolerance>1e-16</error_tolerance-->
</eigen>
<petsc>
<prefix>sd</prefix>
<parameters>-sd_ksp_type cg -sd_pc_type bjacobi -sd_ksp_rtol 1e-16 -sd_ksp_max_it 10000</parameters>
</petsc>
</linear_solver>
</linear_solvers>
<test_definition>
<vtkdiff>
<file>equilibrium_restart_pcs_0_ts_1_t_2.000000.vtu</file>
<field>displacement</field>
<absolute_tolerance>1e-15</absolute_tolerance>
<relative_tolerance>0</relative_tolerance>
</vtkdiff>
<vtkdiff>
<file>equilibrium_restart_pcs_0_ts_1_t_2.000000.vtu</file>
<field>sigma</field>
<absolute_tolerance>1e-7</absolute_tolerance>
<relative_tolerance>0</relative_tolerance>
</vtkdiff>
</test_definition>
</OpenGeoSysProject>
Git LFS file not shown
Loading

0 comments on commit 29ce4b5

Please sign in to comment.