Commit 177209a8 authored by KeitaYoshioka's avatar KeitaYoshioka
Browse files

Another model of phasefield (AT1) has been added.

parent 6a57cea2
Choice for AT1 model (phasefield_model="AT1") or AT2 model (phasefield_model="AT2") of phase-field.
For further information, refer to:
Tann{\'e'} E et al (2018). Crack nucleation in variational phase-field models of brittle fracture. J Mech Phys Solids, 110:80-99
\cite tanne2018crack
......@@ -23,6 +23,16 @@
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 ",
......
......@@ -171,22 +171,37 @@ std::unique_ptr<Process> createPhaseFieldProcess(
const bool hydro_crack = (crack_scheme && (*crack_scheme == "propagating"));
const bool crack_pressure = crack_scheme.has_value();
auto pf_irrv_read =
auto const irreversible_threshold =
//! \ogs_file_param{prj__processes__process__PHASE_FIELD__irreversible_threshold}
config.getConfigParameterOptional<double>("irreversible_threshold");
config.getConfigParameter<double>("irreversible_threshold", 0.05);
double irreversible_threshold = 0.05;
if (pf_irrv_read)
auto const phasefield_model = [&]
{
irreversible_threshold = *pf_irrv_read;
}
auto const phasefield_model_string =
//! \ogs_file_param{prj__processes__process__PHASE_FIELD__phasefield_model}
config.getConfigParameter<std::string>("phasefield_model");
if (phasefield_model_string == "AT1")
{
return PhaseFieldModel::AT1;
}
else if (phasefield_model_string == "AT2")
{
return PhaseFieldModel::AT2;
}
OGS_FATAL(
"phasefield_model must be 'AT1' or 'AT2' but '{:s}' "
"was given",
phasefield_model_string.c_str());
}();
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};
crack_pressure, irreversible_threshold,
phasefield_model};
SecondaryVariableCollection secondary_variables;
......
......@@ -170,25 +170,22 @@ void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
double const gc = _process_data.crack_resistance(t, x_position)[0];
double const ls = _process_data.crack_length_scale(t, x_position)[0];
// for hydro crack, u is rescaled.
if (_process_data.hydro_crack)
{
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;
auto const x_coord =
NumLib::interpolateXCoordinate<ShapeFunction,
ShapeMatricesType>(_element, N);
auto const& B = LinearBMatrix::computeBMatrix<
DisplacementDim, ShapeFunction::NPOINTS,
typename BMatricesType::BMatrixType>(dNdx, N, x_coord,
_is_axially_symmetric);
auto& eps = _ip_data[ip].eps;
eps.noalias() = B * u;
_ip_data[ip].updateConstitutiveRelation(t, x_position, dt, u,
degradation);
}
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;
auto const x_coord =
NumLib::interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
_element, N);
auto const& B =
LinearBMatrix::computeBMatrix<DisplacementDim,
ShapeFunction::NPOINTS,
typename BMatricesType::BMatrixType>(
dNdx, N, x_coord, _is_axially_symmetric);
auto& eps = _ip_data[ip].eps;
eps.noalias() = B * u;
_ip_data[ip].updateConstitutiveRelation(t, x_position, dt, u,
degradation);
auto const& strain_energy_tensile = _ip_data[ip].strain_energy_tensile;
......@@ -209,16 +206,42 @@ void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
}
local_Jac.noalias() +=
(2 * N.transpose() * N * strain_energy_tensile +
gc * (N.transpose() * N / ls + dNdx.transpose() * dNdx * ls)) *
w;
(2 * N.transpose() * N * strain_energy_tensile) * w;
local_rhs.noalias() -=
(N.transpose() * N * d * 2 * strain_energy_tensile +
gc * ((N.transpose() * N / ls + dNdx.transpose() * dNdx * ls) * d -
N.transpose() / ls) -
(N.transpose() * N * d * 2 * strain_energy_tensile -
local_pressure * dNdx.transpose() * N_u * u) *
w;
switch (_process_data.phasefield_model)
{
case PhaseFieldModel::AT1:
{
local_Jac.noalias() +=
gc * 0.75 * dNdx.transpose() * dNdx * ls * w;
local_rhs.noalias() -=
gc *
(-0.375 * N.transpose() / ls +
0.75 * dNdx.transpose() * dNdx * ls * d) *
w;
break;
}
case PhaseFieldModel::AT2:
{
local_Jac.noalias() +=
gc *
(N.transpose() * N / ls + dNdx.transpose() * dNdx * ls) * w;
local_rhs.noalias() -=
gc *
((N.transpose() * N / ls + dNdx.transpose() * dNdx * ls) *
d -
N.transpose() / ls) *
w;
break;
}
}
}
}
......
......@@ -135,6 +135,8 @@ public:
using PhaseFieldMatrix =
typename ShapeMatricesType::template MatrixType<phasefield_size,
phasefield_size>;
using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType;
using NodalVectorType = typename ShapeMatricesType::NodalVectorType;
using IpData =
IntegrationPointData<BMatricesType, ShapeMatricesType, DisplacementDim>;
......
......@@ -32,6 +32,12 @@ struct Parameter;
namespace PhaseField
{
enum class PhaseFieldModel
{
AT1,
AT2
};
template <int DisplacementDim>
struct PhaseFieldProcessData
{
......@@ -47,7 +53,8 @@ struct PhaseFieldProcessData
Eigen::Matrix<double, DisplacementDim, 1> const specific_body_force;
bool hydro_crack = false;
bool crack_pressure = false;
double irreversible_threshold = 0.05;
double irreversible_threshold;
PhaseFieldModel phasefield_model;
double const unity_pressure = 1.0;
double pressure = 0.0;
......
AddTest(
NAME PhaseField_3D_beam_tens_AT1_iso
PATH PhaseField/beam
EXECUTABLE ogs
EXECUTABLE_ARGS AT1_iso_tensile.prj
WRAPPER mpirun
WRAPPER_ARGS -np 1
TESTER vtkdiff
REQUIREMENTS OGS_USE_MPI
RUNTIME 18
DIFF_DATA
expected_AT1_iso_tension_ts_10_t_1_000000_0.vtu AT1_iso_tension_ts_10_t_1_000000_0.vtu displacement displacement 1e-5 0
expected_AT1_iso_tension_ts_10_t_1_000000_0.vtu AT1_iso_tension_ts_10_t_1_000000_0.vtu phasefield phasefield 1e-6 0
)
AddTest(
NAME PhaseField_3D_beam_tens_AT2_iso
PATH PhaseField/beam
......@@ -12,3 +27,4 @@ AddTest(
expected_AT2_iso_tension_ts_10_t_1_000000_0.vtu AT2_iso_tension_ts_10_t_1_000000_0.vtu displacement displacement 1e-5 0
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
)
<?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>
<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_iso_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>
......@@ -12,6 +12,8 @@
<type>PHASE_FIELD</type>
<coupling_scheme>staggered</coupling_scheme>
<integration_order>2</integration_order>
<phasefield_model>AT2</phasefield_model>
<irreversible_threshold>0.05</irreversible_threshold>
<constitutive_relation>
<type>LinearElasticIsotropic</type>
<youngs_modulus>E</youngs_modulus>
......
Markdown is supported
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