List of tutorials
All tutorials are written assuming that you are reasonably familiar with MOOSE. If you find most of the tutorials difficult to follow, please refer to the official MOOSE website for learning resources.
Tutorial 8: Ductile fracture
For simplicity, we will consider a single element 3 dimensional homogeneous cube under tension. For this tutorial we will apply J2 plasticity with a linear hardening model, but in the tutorial folder you will find examples with other configurations.
elastoplasticity.i
: The displacement subproblem
Similar to previous tutorials, we need to define parameters. We define additional parameters related to plasticity.
E = 6.88e4
nu = 0.33
K = '${fparse E/3/(1-2*nu)}'
G = '${fparse E/2/(1+nu)}'
sigma_y = 320
H = 3.44e3
Where H
is hardening modulus and sigma_y
is the yield stress.
GeneratedMeshGenerator
is utilized to create a 3 dimensional cube.
[Mesh]
[gmg]
type = GeneratedMeshGenerator
dim = 3
xmax = ${a}
ymax = ${a}
zmax = ${a}
[]
[]
We are now transferring an additional variable to the fracture multi-app, psip_active. This is the active plastic energy. In this model active plastic energy contributes to the fracture driving energy.
[Transfers]
[from_d]
type = MultiAppCopyTransfer
from_multi_app = 'fracture'
variable = d
source_variable = d
[]
[to_psie_active]
type = MultiAppCopyTransfer
to_multi_app = 'fracture'
variable = psie_active
source_variable = psie_active
[]
[to_psip_active]
type = MultiAppCopyTransfer
to_multi_app = 'fracture'
variable = psip_active
source_variable = psip_active
[]
[]
We have an additional auxiliary variable defined, F
. This is a consequence of utilizing large deformation models. Previous tutorials utilized small deformation models.
[AuxVariables]
[d]
[]
[stress]
order = CONSTANT
family = MONOMIAL
[]
[F]
order = CONSTANT
family = MONOMIAL
[]
[]
Plasticity is added through constitutive updates in the Materials block.
[Materials]
[bulk_properties]
type = ADGenericConstantMaterial
prop_names = 'K G l Gc psic sigma_y H'
prop_values = '${K} ${G} ${l} ${Gc} ${psic} ${sigma_y} ${H}'
[]
[degradation]
type = RationalDegradationFunction
property_name = g
expression = (1-d)^p/((1-d)^p+(Gc/psic*xi/c0/l)*d*(1+a2*d+a2*a3*d^2))*(1-eta)+eta
phase_field = d
material_property_names = 'Gc psic xi c0 l '
parameter_names = 'p a2 a3 eta '
parameter_values = '2 -0.5 0 1e-6'
[]
[crack_geometric]
type = CrackGeometricFunction
property_name = alpha
expression = 'd'
phase_field = d
[]
[defgrad]
type = ComputeDeformationGradient
[]
[hencky]
type = HenckyIsotropicElasticity
bulk_modulus = K
shear_modulus = G
phase_field = d
degradation_function = g
output_properties = 'psie_active'
outputs = exodus
[]
[linear_hardening]
type = LinearHardening
phase_field = d
yield_stress = sigma_y
hardening_modulus = H
degradation_function = g
output_properties = 'psip_active'
outputs = exodus
[]
[J2]
type = LargeDeformationJ2Plasticity
hardening_model = linear_hardening
[]
[stress]
type = ComputeLargeDeformationStress
elasticity_model = hencky
plasticity_model = J2
[]
[]
Here we add the blocks for HenckyIsotropicElasticity
, linearhardening
, LargeDeformationJ2Plasticity
, and ComputeLargeDeformationStress
.
LargeDeformationJ2Plasticity
: This utilizes the selected hardening model to calculate the value of the yield surface and determine whether the material has begun plastic deformation.
HenckyIsotropicElasticity:
This calculates elastic stress using the hencky model. Returns the Mandel stress to ComputeLargeDeformationStress
LinearHardening:
This performs plastic hardening following a linear isotropic law. It returns plastic energy to LargeDeformationJ2Plasticity
ComputeLargeDeformationStress:
This computes stress given an elasticity model and a plasticity model. It uses the elasticity model to compute Cauchy stress and the plasticity model to update the plastic deformation to bring the stress state back unto the yield surface.
Modifications
The elasticity model and hardening law used can be changed.
Currently available hyperelasticity models
HenckyIsotropicElasticity
CNHIsotropicElasticity
Currently available hardening models
ArrheniusLawHardening
LinearHardening
PowerLawHardening
fracture.i
: The phase-field subproblem
There are two minor changes to be made to the fracture input file seen in previous tutorials.
First, we add an auxiliary variable for active plastic energy. This is to receive the value from the main app.
[AuxVariables]
[bounds_dummy]
[]
[psie_active]
order = CONSTANT
family = MONOMIAL
[]
[psip_active]
order = CONSTANT
family = MONOMIAL
[]
[]
Second, we add the active plastic energy into the fracture driving energy
[Materials]
[psi]
type = ADDerivativeParsedMaterial
property_name = psi
expression = 'alpha*Gc/c0/l+g*(psie_active+psip_active)'
coupled_variables = 'd psie_active psip_active'
material_property_names = 'alpha(d) g(d) Gc c0 l'
derivative_order = 1
[]
[]
Alternating Minimization
Like in previous tutorials, we rely on on MOOSE's Picard iteration to create an alternating minimization scheme.
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options_iname = '-pc_type -snes_type'
petsc_options_value = 'lu vinewtonrsls'
nl_rel_tol = 1e-08
nl_abs_tol = 1e-10
nl_max_its = 50
automatic_scaling = true
abort_on_solve_fail = true
[]