Skip to content
Snippets Groups Projects
Commit fa5c0498 authored by Dmitry Yu. Naumov's avatar Dmitry Yu. Naumov
Browse files

Merge branch 'PhaseField' into 'master' 

Additional material model for Phase field

See merge request ogs/ogs!3850
parents a9786e21 b6863e4f
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 510 additions and 54 deletions
Documentation/ProjectFile/prj/processes/process/PHASE_FIELD/t_energy_split_model.md 0 → 100644
+ 3
− 0
View file @ fa5c0498
Implemented energy split models are:
- isotropic model \cite bourdin2000numerical, (energy_split_model="Isotropic")
- volumetric-deviatoric split model \cite amor2009regularized (energy_split_model="VolumetricDeviatoric")
Documentation/bibliography/ogs.bib
+ 0
− 10
View file @ fa5c0498
......@@ -23,16 +23,6 @@
category={Methods},
}
@article{tanne2018crack,
title={Crack nucleation in variational phase-field models of brittle fracture},
author={Tann{\'e}, Erwan and Li, Tianyi and Bourdin, Blaise and Marigo, J-J and Maurini, Corrado},
journal={Journal of the Mechanics and Physics of Solids},
volume={110},
pages={80--99},
year={2018},
publisher={Elsevier}
}
@article{Parisio2017,
title = "Plastic-damage modeling of saturated quasi-brittle shales ",
journal = "International Journal of Rock Mechanics and Mining Sciences ",
......
Documentation/bibliography/other.bib
+ 33
− 0
View file @ fa5c0498
......@@ -32,6 +32,39 @@ author = {J. Buchwald and S. Kaiser and O. Kolditz and T. Nagel}
Publisher={Springer}
}
@article{tanne2018crack,
title={Crack nucleation in variational phase-field models of brittle fracture},
author={Tann{\'e}, Erwan and Li, Tianyi and Bourdin, Blaise and Marigo, J-J and Maurini, Corrado},
journal={Journal of the Mechanics and Physics of Solids},
volume={110},
pages={80--99},
year={2018},
publisher={Elsevier}
}
@article{bourdin2000numerical,
title={Numerical experiments in revisited brittle fracture},
author={Bourdin, Blaise and Francfort, Gilles A and Marigo, Jean-Jacques},
journal={Journal of the Mechanics and Physics of Solids},
volume={48},
number={4},
pages={797--826},
year={2000},
publisher={Elsevier}
}
@article{amor2009regularized,
title={Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments},
author={Amor, Hanen and Marigo, Jean-Jacques and Maurini, Corrado},
journal={Journal of the Mechanics and Physics of Solids},
volume={57},
number={8},
pages={1209--1229},
year={2009},
publisher={Elsevier}
}
@Book{Ericson:2004:RCD:1121584,
Title = {Real-Time Collision Detection (The Morgan Kaufmann Series in Interactive 3-D Technology)},
Author = {Ericson, Christer},
......
MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.h
+ 52
− 3
View file @ fa5c0498
......@@ -32,10 +32,9 @@ std::tuple<MathLib::KelvinVector::KelvinVectorType<
DisplacementDim> /* C_tensile */,
MathLib::KelvinVector::KelvinMatrixType<
DisplacementDim> /* C_compressive */,
double /* strain_energy_tensile */,
double /* elastic_energy */
double /* strain_energy_tensile */, double /* elastic_energy */
>
calculateDegradedStress(
calculateVolDevDegradedStress(
double const degradation,
double const bulk_modulus,
double const mu,
......@@ -89,6 +88,56 @@ calculateDegradedStress(
return std::make_tuple(sigma_real, sigma_tensile, C_tensile, C_compressive,
strain_energy_tensile, elastic_energy);
}
template <int DisplacementDim>
std::tuple<MathLib::KelvinVector::KelvinVectorType<
DisplacementDim> /* sigma_real */,
MathLib::KelvinVector::KelvinVectorType<
DisplacementDim> /* sigma_tensile */,
MathLib::KelvinVector::KelvinMatrixType<
DisplacementDim> /* C_tensile */,
MathLib::KelvinVector::KelvinMatrixType<
DisplacementDim> /* C_compressive */,
double /* strain_energy_tensile */, double /* elastic_energy */
>
calculateIsotropicDegradedStress(
double const degradation,
double const bulk_modulus,
double const mu,
MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const& eps)
{
static int const KelvinVectorSize =
MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
using KelvinVector =
MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
using KelvinMatrix =
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>;
using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
// calculation of deviatoric parts
auto const& P_dev = Invariants::deviatoric_projection;
KelvinVector const epsd_curr = P_dev * eps;
// Hydrostatic part for the stress and the tangent.
double const eps_curr_trace = Invariants::trace(eps);
KelvinMatrix C_tensile = KelvinMatrix::Zero();
KelvinMatrix C_compressive = KelvinMatrix::Zero();
double const strain_energy_tensile =
bulk_modulus / 2 * eps_curr_trace * eps_curr_trace +
mu * epsd_curr.transpose() * epsd_curr;
double const elastic_energy = degradation * strain_energy_tensile;
KelvinVector const sigma_tensile =
bulk_modulus * eps_curr_trace * Invariants::identity2 +
2 * mu * epsd_curr;
C_tensile.template topLeftCorner<3, 3>().setConstant(bulk_modulus);
C_tensile.noalias() += 2 * mu * P_dev * KelvinMatrix::Identity();
KelvinVector const sigma_real = degradation * sigma_tensile;
return std::make_tuple(sigma_real, sigma_tensile, C_tensile, C_compressive,
strain_energy_tensile, elastic_energy);
}
} // namespace Phasefield
} // namespace Solids
} // namespace MaterialLib
ProcessLib/PhaseField/CreatePhaseFieldProcess.cpp
+ 21
− 1
View file @ fa5c0498
......@@ -195,13 +195,33 @@ std::unique_ptr<Process> createPhaseFieldProcess(
phasefield_model_string.c_str());
}();
auto const energy_split_model = [&]
{
auto const energy_split_model_string =
//! \ogs_file_param{prj__processes__process__PHASE_FIELD__energy_split_model}
config.getConfigParameter<std::string>("energy_split_model");
if (energy_split_model_string == "Isotropic")
{
return EnergySplitModel::Isotropic;
}
else if (energy_split_model_string == "VolumetricDeviatoric")
{
return EnergySplitModel::VolDev;
}
OGS_FATAL(
"energy_split_model must be 'Isotropic' or 'VolumetricDeviatoric' "
"but '{:s}' was given",
energy_split_model_string);
}();
PhaseFieldProcessData<DisplacementDim> process_data{
materialIDs(mesh), std::move(solid_constitutive_relations),
residual_stiffness, crack_resistance,
crack_length_scale, solid_density,
specific_body_force, hydro_crack,
crack_pressure, irreversible_threshold,
phasefield_model};
phasefield_model, energy_split_model};
SecondaryVariableCollection secondary_variables;
......
ProcessLib/PhaseField/PhaseFieldFEM-impl.h
+ 6
− 4
View file @ fa5c0498
......@@ -95,8 +95,9 @@ void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
double const k = _process_data.residual_stiffness(t, x_position)[0];
double const d_ip = N.dot(d);
double const degradation = d_ip * d_ip * (1 - k) + k;
_ip_data[ip].updateConstitutiveRelation(t, x_position, dt, u,
degradation);
_ip_data[ip].updateConstitutiveRelation(
t, x_position, dt, u, degradation,
_process_data.energy_split_model);
auto& sigma = _ip_data[ip].sigma;
auto const& C_tensile = _ip_data[ip].C_tensile;
......@@ -184,8 +185,9 @@ void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
auto& eps = _ip_data[ip].eps;
eps.noalias() = B * u;
_ip_data[ip].updateConstitutiveRelation(t, x_position, dt, u,
degradation);
_ip_data[ip].updateConstitutiveRelation(
t, x_position, dt, u, degradation,
_process_data.energy_split_model);
auto const& strain_energy_tensile = _ip_data[ip].strain_energy_tensile;
......
ProcessLib/PhaseField/PhaseFieldFEM.h
+ 27
− 12
View file @ fa5c0498
......@@ -74,7 +74,8 @@ struct IntegrationPointData final
ParameterLib::SpatialPosition const& x,
double const /*dt*/,
DisplacementVectorType const& /*u*/,
double const degradation)
double const degradation,
EnergySplitModel const energy_split_model)
{
auto linear_elastic_mp =
static_cast<MaterialLib::Solids::LinearElasticIsotropic<
......@@ -84,22 +85,35 @@ struct IntegrationPointData final
auto const bulk_modulus = linear_elastic_mp.bulk_modulus(t, x);
auto const mu = linear_elastic_mp.mu(t, x);
std::tie(sigma, sigma_tensile, C_tensile, C_compressive,
strain_energy_tensile, elastic_energy) =
MaterialLib::Solids::Phasefield::calculateDegradedStress<
DisplacementDim>(degradation, bulk_modulus, mu, eps);
switch (energy_split_model)
{
case EnergySplitModel::Isotropic:
{
std::tie(sigma, sigma_tensile, C_tensile, C_compressive,
strain_energy_tensile, elastic_energy) =
MaterialLib::Solids::Phasefield::
calculateIsotropicDegradedStress<DisplacementDim>(
degradation, bulk_modulus, mu, eps);
break;
}
case EnergySplitModel::VolDev:
{
std::tie(sigma, sigma_tensile, C_tensile, C_compressive,
strain_energy_tensile, elastic_energy) =
MaterialLib::Solids::Phasefield::
calculateVolDevDegradedStress<DisplacementDim>(
degradation, bulk_modulus, mu, eps);
break;
}
}
history_variable =
std::max(history_variable_prev, strain_energy_tensile);
}
EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
static constexpr int kelvin_vector_size =
MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
};
/// Used for the extrapolation of the integration point values. It is ordered
/// (and stored) by integration points.
/// Used for the extrapolation of the integration point values. It is
/// ordered (and stored) by integration points.
template <typename ShapeMatrixType>
struct SecondaryData
{
......@@ -242,7 +256,8 @@ public:
}
void postTimestepConcrete(Eigen::VectorXd const& /*local_x*/,
double const /*t*/, double const /*dt*/) override
double const /*t*/,
double const /*dt*/) override
{
unsigned const n_integration_points =
_integration_method.getNumberOfPoints();
......
ProcessLib/PhaseField/PhaseFieldProcessData.h
+ 7
− 0
View file @ fa5c0498
......@@ -38,6 +38,12 @@ enum class PhaseFieldModel
AT2
};
enum class EnergySplitModel
{
Isotropic,
VolDev
};
template <int DisplacementDim>
struct PhaseFieldProcessData
{
......@@ -55,6 +61,7 @@ struct PhaseFieldProcessData
bool crack_pressure = false;
double irreversible_threshold;
PhaseFieldModel phasefield_model;
EnergySplitModel energy_split_model;
double const unity_pressure = 1.0;
double pressure = 0.0;
......
ProcessLib/PhaseField/Tests.cmake
+ 14
− 0
View file @ fa5c0498
......@@ -28,3 +28,17 @@ AddTest(
expected_AT2_iso_tension_ts_10_t_1_000000_0.vtu AT2_iso_tension_ts_10_t_1_000000_0.vtu phasefield phasefield 1e-6 0
)
AddTest(
NAME PhaseField_3D_beam_tens_AT1_vd
PATH PhaseField/beam
EXECUTABLE ogs
EXECUTABLE_ARGS AT1_vd_tensile.prj
WRAPPER mpirun
WRAPPER_ARGS -np 1
TESTER vtkdiff
REQUIREMENTS OGS_USE_MPI
RUNTIME 18
DIFF_DATA
expected_AT1_vd_tension_ts_10_t_1_000000_0.vtu AT1_vd_tension_ts_10_t_1_000000_0.vtu displacement displacement 1e-5 0
expected_AT1_vd_tension_ts_10_t_1_000000_0.vtu AT1_vd_tension_ts_10_t_1_000000_0.vtu phasefield phasefield 1e-6 0
)
ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM.h
+ 1
− 1
View file @ fa5c0498
......@@ -92,7 +92,7 @@ struct IntegrationPointData final
std::tie(sigma, sigma_tensile, C_tensile, C_compressive,
strain_energy_tensile, elastic_energy) =
MaterialLib::Solids::Phasefield::calculateDegradedStress<
MaterialLib::Solids::Phasefield::calculateVolDevDegradedStress<
DisplacementDim>(degradation, bulk_modulus, mu, eps_m);
}
EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
......
Tests/Data/PhaseField/beam/AT1_iso_tensile.prj
+ 1
− 0
View file @ fa5c0498
......@@ -13,6 +13,7 @@
<coupling_scheme>staggered</coupling_scheme>
<integration_order>2</integration_order>
<phasefield_model>AT1</phasefield_model>
<energy_split_model>Isotropic</energy_split_model>
<irreversible_threshold>0.05</irreversible_threshold>
<constitutive_relation>
<type>LinearElasticIsotropic</type>
......
Tests/Data/PhaseField/beam/AT1_vd_tensile.prj 0 → 100644
+ 290
− 0
View file @ fa5c0498
<?xml version="1.0" encoding="ISO-8859-1"?>
<OpenGeoSysProject>
<!-- element size == 0.01 -->
<meshes>
<mesh>bar.vtu</mesh>
<mesh>bar_left.vtu</mesh>
<mesh>bar_right.vtu</mesh>
</meshes>
<processes>
<process>
<name>PhaseField</name>
<type>PHASE_FIELD</type>
<coupling_scheme>staggered</coupling_scheme>
<integration_order>2</integration_order>
<phasefield_model>AT1</phasefield_model>
<energy_split_model>VolumetricDeviatoric</energy_split_model>
<irreversible_threshold>0.05</irreversible_threshold>
<constitutive_relation>
<type>LinearElasticIsotropic</type>
<youngs_modulus>E</youngs_modulus>
<poissons_ratio>nu</poissons_ratio>
</constitutive_relation>
<phasefield_parameters>
<residual_stiffness>k</residual_stiffness>
<crack_resistance>gc</crack_resistance>
<crack_length_scale>ls</crack_length_scale>
</phasefield_parameters>
<solid_density>rho_sr</solid_density>
<process_variables>
<phasefield>phasefield</phasefield>
<displacement>displacement</displacement>
</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>
<specific_body_force>0 0 0</specific_body_force>
</process>
</processes>
<time_loop>
<global_process_coupling>
<max_iter> 1000 </max_iter>
<convergence_criteria>
<!-- convergence criterion for the first process -->
<convergence_criterion>
<type>DeltaX</type>
<norm_type>INFINITY_N</norm_type>
<abstol>1.e-1</abstol>
<reltol>1.e-1</reltol>
</convergence_criterion>
<!-- convergence criterion for the second process -->
<convergence_criterion>
<type>DeltaX</type>
<norm_type>INFINITY_N</norm_type>
<abstol>1.e-4</abstol>
<reltol>1.e-10</reltol>
</convergence_criterion>
</convergence_criteria>
</global_process_coupling>
<processes>
<process ref="PhaseField">
<nonlinear_solver>basic_newton_u</nonlinear_solver>
<convergence_criterion>
<type>DeltaX</type>
<norm_type>NORM2</norm_type>
<reltol>1.e-14</reltol>
</convergence_criterion>
<time_discretization>
<type>BackwardEuler</type>
</time_discretization>
<time_stepping>
<type>FixedTimeStepping</type>
<t_initial>0</t_initial>
<t_end>1.0</t_end>
<timesteps>
<pair>
<repeat>1</repeat>
<delta_t>1.e-1</delta_t>
</pair>
</timesteps>
</time_stepping>
</process>
<process ref="PhaseField">
<nonlinear_solver>petsc_snes</nonlinear_solver>
<convergence_criterion>
<type>DeltaX</type>
<norm_type>NORM2</norm_type>
<reltol>1.e-14</reltol>
</convergence_criterion>
<time_discretization>
<type>BackwardEuler</type>
</time_discretization>
<time_stepping>
<type>FixedTimeStepping</type>
<t_initial>0</t_initial>
<t_end>1.0</t_end>
<timesteps>
<pair>
<repeat>1</repeat>
<delta_t>1.e-1</delta_t>
</pair>
</timesteps>
</time_stepping>
</process>
</processes>
<output>
<variables>
<variable>displacement</variable>
<variable>phasefield</variable>
<variable>sigma</variable>
<variable>epsilon</variable>
</variables>
<type>VTK</type>
<prefix>AT1_vd_tension</prefix>
<timesteps>
<pair>
<repeat>10</repeat>
<each_steps>1</each_steps>
</pair>
<pair>
<repeat>10</repeat>
<each_steps>1</each_steps>
</pair>
</timesteps>
</output>
</time_loop>
<parameters>
<!-- Mechanics -->
<parameter>
<name>E</name>
<type>Constant</type>
<value>1.0</value>
</parameter>
<parameter>
<name>nu</name>
<type>Constant</type>
<value>0.15</value>
</parameter>
<parameter>
<name>k</name>
<type>Constant</type>
<value>1e-8</value>
</parameter>
<parameter>
<name>gc</name>
<type>Constant</type>
<value>1.0</value>
</parameter>
<parameter>
<name>ls</name>
<type>Constant</type>
<value>0.02</value>
</parameter>
<parameter>
<name>rho_sr</name>
<type>Constant</type>
<value>0.0</value>
</parameter>
<parameter>
<name>displacement0</name>
<type>Constant</type>
<values>0 0 0</values>
</parameter>
<parameter>
<name>phasefield_ic</name>
<type>Constant</type>
<value>1</value>
</parameter>
<parameter>
<name>phasefield_bc</name>
<type>Constant</type>
<value>1</value>
</parameter>
<parameter>
<name>dirichlet0</name>
<type>Constant</type>
<value>0</value>
</parameter>
<parameter>
<name>Dirichlet_spatial</name>
<type>Constant</type>
<value>1</value>
</parameter>
<parameter>
<name>dirichlet_right_time</name>
<type>CurveScaled</type>
<curve>dirichlet_time</curve>
<parameter>dirichlet_right</parameter>
</parameter>
<parameter>
<name>dirichlet_right</name>
<type>Constant</type>
<value>5.0</value>
</parameter>
</parameters>
<curves>
<curve>
<name>dirichlet_time</name>
<coords>0 1.0</coords>
<values>0 0.9</values>
</curve>
</curves>
<process_variables>
<process_variable>
<name>phasefield</name>
<components>1</components>
<order>1</order>
<initial_condition>phasefield_ic</initial_condition>
<boundary_conditions>
<boundary_condition>
<geometrical_set>bar</geometrical_set>
<geometry>left</geometry>
<type>Dirichlet</type>
<component>0</component>
<parameter>phasefield_bc</parameter>
</boundary_condition>
<boundary_condition>
<geometrical_set>bar</geometrical_set>
<geometry>right</geometry>
<type>Dirichlet</type>
<component>0</component>
<parameter>phasefield_bc</parameter>
</boundary_condition>
<boundary_condition>
<geometrical_set>bar</geometrical_set>
<geometry>right</geometry>
<type>PhaseFieldIrreversibleDamageOracleBoundaryCondition</type>
<component>0</component>
</boundary_condition>
</boundary_conditions>
</process_variable>
<process_variable>
<name>displacement</name>
<components>3</components>
<order>1</order>
<initial_condition>displacement0</initial_condition>
<boundary_conditions>
<boundary_condition>
<geometrical_set>bar</geometrical_set>
<geometry>left</geometry>
<type>Dirichlet</type>
<component>0</component>
<parameter>dirichlet0</parameter>
</boundary_condition>
<boundary_condition>
<geometrical_set>bar</geometrical_set>
<geometry>right</geometry>
<type>Dirichlet</type>
<component>0</component>
<parameter>dirichlet_right_time</parameter>
</boundary_condition>
</boundary_conditions>
</process_variable>
</process_variables>
<nonlinear_solvers>
<nonlinear_solver>
<name>petsc_snes</name>
<type>PETScSNES</type>
<!-- <type>Newton</type> -->
<max_iter>50</max_iter>
<linear_solver>linear_solver_d</linear_solver>
</nonlinear_solver>
<nonlinear_solver>
<name>basic_newton_u</name>
<type>Newton</type>
<max_iter>200</max_iter>
<linear_solver>linear_solver_u</linear_solver>
</nonlinear_solver>
</nonlinear_solvers>
<linear_solvers>
<linear_solver>
<name>linear_solver_d</name>
<eigen>
<solver_type>BiCGSTAB</solver_type>
<precon_type>ILUT</precon_type>
<max_iteration_step>10000</max_iteration_step>
<error_tolerance>1e-16</error_tolerance>
</eigen>
<petsc>
<parameters>-ksp_type cg -pc_type bjacobi -ksp_atol 1e-10 -ksp_rtol 1e-10 -snes_type vinewtonrsls -snes_linesearch_type l2 -snes_atol 1.e-8 -snes_rtol 1.e-8 -snes_max_it 1000 -snes_monitor </parameters>
</petsc>
</linear_solver>
<linear_solver>
<name>linear_solver_u</name>
<petsc>
<parameters>-ksp_type cg -pc_type bjacobi -ksp_atol 1e-14 -ksp_rtol 1e-14 </parameters>
</petsc>
</linear_solver>
</linear_solvers>
</OpenGeoSysProject>
Tests/Data/PhaseField/beam/AT2_iso_tensile.prj
+ 1
− 0
View file @ fa5c0498
......@@ -13,6 +13,7 @@
<coupling_scheme>staggered</coupling_scheme>
<integration_order>2</integration_order>
<phasefield_model>AT2</phasefield_model>
<energy_split_model>Isotropic</energy_split_model>
<irreversible_threshold>0.05</irreversible_threshold>
<constitutive_relation>
<type>LinearElasticIsotropic</type>
......
Tests/Data/PhaseField/beam/expected_AT1_vd_tension_ts_10_t_1_000000_0.vtu 0 → 100644
+ 31
− 0
View file @ fa5c0498
This diff is collapsed.
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment