Bilinear mixed mode traction separation law

This class implements the bilinear mixed mode traction separation law described in Camanho and Dávila (2002).

Softening onset prediction

The initiation of the softening process is predicted using the quadratic failure criterion given below,

$var element = document.getElementById("moose-equation-682df66f-9297-4b4f-91d9-f40fc8412038");katex.render("\\left(\\frac{\\langle\\tau_3\\rangle}{N}\\right)^2 + \\left(\\frac{\\tau_2}{S} \\right)^2 + \\left(\\frac{\\tau_1}{T} \\right)^2 = 1", element, {displayMode:true,throwOnError:false});$

The total mixed-mode relative displacement $var element = document.getElementById("moose-equation-a0b3d55a-b3a1-43e7-b2d7-b56872374e17");katex.render("\\delta_m", element, {displayMode:false,throwOnError:false});$ is defined as $var element = document.getElementById("moose-equation-1a2f9d63-7ae7-41f9-a0b7-6d568c7f187e");katex.render("\\delta_m = \\sqrt{\\delta_1^2+\\delta_2^2+\\langle\\delta_3\\rangle^2} = \\sqrt{\\delta_{shear}^2+\\langle\\delta_3\\rangle^2}", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-21df1b6c-ea6d-45ae-b087-74f66004d7ca");katex.render("\\delta_{shear}", element, {displayMode:false,throwOnError:false});$ represents the norm of the vector defining the tangential relative displacements of the element.

Using the same penalty stiffness in Modes I, II and III, the tractions before softening onset are: $var element = document.getElementById("moose-equation-914e237f-8528-4a98-a997-b3f887e00a59");katex.render("\\tau_i = K\\delta_i, ~~~\\text{i=1,2,3}", element, {displayMode:true,throwOnError:false});$

Assuming $var element = document.getElementById("moose-equation-39ced233-c7b4-499a-841b-e79a7b01f6a9");katex.render("S=T", element, {displayMode:false,throwOnError:false});$ , the single mode relative displacements at softening onset are: $var element = document.getElementById("moose-equation-779ca502-9617-414d-82c9-e8dc7c848b7b");katex.render("\\delta_3^0=\\frac{N}{K}", element, {displayMode:true,throwOnError:false});$

$var element = document.getElementById("moose-equation-fb509618-1875-4da9-8602-d053b99ea96e");katex.render("\\delta_1^0=\\delta_2^0=\\delta_{shear}^0=\\frac{S}{K}", element, {displayMode:true,throwOnError:false});$

For an opening displacement $var element = document.getElementById("moose-equation-139dff62-059e-48e6-acc0-6aba1f6538ed");katex.render("\\delta_3", element, {displayMode:false,throwOnError:false});$ greater than zero, the mode mixity ratio $var element = document.getElementById("moose-equation-d0671ced-bac8-4007-9baf-cccb1fe3a4a8");katex.render("\\beta", element, {displayMode:false,throwOnError:false});$ is defined as: $var element = document.getElementById("moose-equation-f9c024ec-1d28-42d0-ba18-aab573dfffda");katex.render("\\beta = \\frac{\\delta_{shear}}{\\delta_3}", element, {displayMode:true,throwOnError:false});$

The mixed-mode relative displacement corresponding to the onset of softening $var element = document.getElementById("moose-equation-7dd0a5c1-d61e-4bf6-9618-4d9484f88714");katex.render("\\delta_m^0", element, {displayMode:false,throwOnError:false});$ is given as $var element = document.getElementById("moose-equation-ba3a2615-c8f1-4697-8a44-68a30c39937a");katex.render("\\delta_m^0= \\begin{cases} \\delta_3^0\\delta_1^0\\sqrt{\\frac{1+\\beta^2}{(\\delta_1^0)^2+(\\beta\\delta_3^0)^2}} , & \\delta_3> 0\\\\ \\delta_{shear}^0, & \\delta_3\\leq 0 \\end{cases}", element, {displayMode:true,throwOnError:false});$

Delamination propagation prediction

Power law criterion

The power law criterion is given as $var element = document.getElementById("moose-equation-46a00aac-7b00-447c-970d-6faad6d8b997");katex.render("\\left(\\frac{G_{I}}{G_{IC}}\\right)^{\\alpha} + \\left(\\frac{G_{II}}{G_{IIC}}\\right)^{\\alpha} = 1", element, {displayMode:true,throwOnError:false});$

The mixed-mode displacements corresponding to total decohesion is given as: $var element = document.getElementById("moose-equation-7c742873-7588-4399-9e01-8b288b9b378a");katex.render("\\delta_m^f= \\begin{cases} \\ \\frac{2(1+\\beta^2)}{K\\delta_m^0}\\left[\\left(\\frac{1}{G_{IC}}\\right)^{\\alpha} + \\left(\\frac{\\beta^2}{G_{IIC}}\\right)^{\\alpha}\\right]^{-1/\\alpha} , & \\delta_3> 0\\\\ \\ \\sqrt{(\\delta_1^f)^2+(\\delta_2^f)^2}, & \\delta_3\\leq 0 \\end{cases}", element, {displayMode:true,throwOnError:false});$

B-K criterion

The mixed-mode criterion proposed by Benzeggagh and Kenane is given as (B-K criterion): $var element = document.getElementById("moose-equation-a7851ce2-3a7c-44fb-b09d-b57ff8afee74");katex.render("G_{IC}+(G_{IIC}-G_{IC})\\left(\\frac{G_{shear}}{G_T}\\right)^{\\eta} = G_C~\\text{with}~G_T=G_I+G_{shear}", element, {displayMode:true,throwOnError:false});$

The mixed-mode displacements corresponding to total decohesion is given as: $var element = document.getElementById("moose-equation-fefb0bfe-a28c-4392-8f1c-ad3284d6a166");katex.render("\\delta_m^f= \\begin{cases} \\ \\frac{2}{K\\delta_m^0}\\left[G_{IC}+(G_{IIC}-G_{IC}\\left(\\frac{\\beta^2}{1+\\beta^2}\\right)^\\eta)\\right] , & \\delta_3> 0\\\\ \\ \\sqrt{(\\delta_1^f)^2+(\\delta_2^f)^2}, & \\delta_3\\leq 0 \\end{cases}", element, {displayMode:true,throwOnError:false});$

Constitutive equation for mixed-mode loading

The constitutive equation for mixed-mode loading is given as $var element = document.getElementById("moose-equation-0a8540ff-5570-400f-a14e-8a793a35afc6");katex.render("\\tau_s = D_{sr}\\delta_r", element, {displayMode:true,throwOnError:false});$

$(1)var element = document.getElementById("moose-equation-f2ebcb49-9606-427d-895a-0909a706bd62");katex.render("D_{sr} = \\begin{cases} \\bar{\\delta}_{sr}K, \\text{if}~\\delta_m^{max} \\leq \\delta_m^0 \\\\ \\bar{\\delta}_{sr}\\left[(1-d)K+Kd\\bar{\\delta}_{s3}\\frac{\\langle-\\delta_3\\rangle}{-\\delta_3}\\right], \\text{if}~\\delta_m^0 < \\delta_m^{max} \\leq \\delta_m^f \\\\ \\bar{\\delta}_{s3}\\bar{\\delta}_{3r}\\frac{\\langle-\\delta_3\\rangle}{-\\delta_3}K, \\text{if}~\\delta_m^{max} \\geq \\delta_m^f \\end{cases} ", element, {displayMode:true,throwOnError:false});$

$(2)var element = document.getElementById("moose-equation-a6c7f985-c540-4550-ac35-28a8c29d55fd");katex.render("d= \\frac{\\delta_m^f(\\delta_m^{max}-\\delta_m^0)}{\\delta_m^{max}(\\delta_m^f-\\delta_m^0)}, d\\in[0,1] ", element, {displayMode:true,throwOnError:false});$

Solver options

Viscous regularization

Cohesive zone models exhibiting softening behavior and stiffness degradation often lead to convergence difficulties in an implicit solver. The traction-separation laws can be regularized using viscosity. The viscous damage variable $var element = document.getElementById("moose-equation-04053892-a807-4041-8d78-21554e7ae585");katex.render("d_v", element, {displayMode:false,throwOnError:false});$ is defined by $var element = document.getElementById("moose-equation-8138639e-2bd7-49ef-a7a0-909d9a7ce6b4");katex.render("\\dot{d_v}=\\frac{1}{\\mu}(d-d_v)", element, {displayMode:true,throwOnError:false});$ where $var element = document.getElementById("moose-equation-e7e22e7d-847c-43a2-af1a-135eecb8fde5");katex.render("\\mu", element, {displayMode:false,throwOnError:false});$ is the viscosity parameter representing the relaxation time of the viscous system. An analytical expression of $var element = document.getElementById("moose-equation-63a14c68-086d-4318-bd40-d1f2fc31e536");katex.render("d_v", element, {displayMode:false,throwOnError:false});$ can be obtained by using the backward Euler method. With viscous regularization, the $var element = document.getElementById("moose-equation-88130bd0-a5fa-444e-b421-dbd7e933220d");katex.render("d", element, {displayMode:false,throwOnError:false});$ will be replaced by $var element = document.getElementById("moose-equation-05344fea-a472-40cd-a3c9-9f12f2e7179c");katex.render("d_v", element, {displayMode:false,throwOnError:false});$ in Eq. (1) to compute traction.

Lag separation state

It is typically useful to improve convergence by lagging the separation state. When lag_seperation_state = true, the $var element = document.getElementById("moose-equation-ae8e23b5-caa1-4364-8875-f86fda4a2fc9");katex.render("\\delta_3", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-ac968cd0-24ae-4c4c-b5c2-89109d829a31");katex.render("\\delta_m^{max}", element, {displayMode:false,throwOnError:false});$ , $var element = document.getElementById("moose-equation-31dc9fff-d78e-465a-917f-b480460cff5f");katex.render("\\delta_m^0", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-0982a3f0-9352-4340-a82a-797b5340a2f7");katex.render("\\delta_m^f", element, {displayMode:false,throwOnError:false});$ will be replaced by their old values from previous time step.

Use Regularized Heavyside Function

The step (heavyside) function $var element = document.getElementById("moose-equation-9c5965c0-cb09-4993-ba24-3e434f9ee0ce");katex.render("\\frac{\\langle-\\delta_3\\rangle}{-\\delta_3}", element, {displayMode:false,throwOnError:false});$ in Eq. (1) usually makes convergence bad. In the code, we replaced it with the regularized heavside function which provides a C0 continuity. The regularization parameter can be set by alpha parameter.

Examples

[czm]
  type = BiLinearMixedModeTraction
  boundary = 'interface'
  penalty_stiffness = 1e6
  GI_c = 1e3
  GII_c = 1e2
  normal_strength = 1e4
  shear_strength = 1e3
  displacements = 'disp_x disp_y'
  eta = 2.2
  viscosity = 1e-3
[]

References

Pedro P. Camanho and Carlos G. Dávila. Mixed-mode decohesion finite elements for the simulation of delamination in composite materials. Technical Report NASA/TM-2002-211737, National Aeronautics and Space Administration, 2002.

@techreport{Camanho2002,
    author = "Camanho, Pedro P. and D{\'a}vila, Carlos G.",
    title = "Mixed-Mode Decohesion Finite Elements for the Simulation of Delamination in Composite Materials",
    year = "2002",
    number = "NASA/TM-2002-211737",
    institution = "National Aeronautics and Space Administration"
}

(moose/modules/tensor_mechanics/test/tests/cohesive_zone_model/bilinear_mixed.i)

[Mesh]
  [msh]
    type = GeneratedMeshGenerator
    dim = 2
    xmax = 1
    ymax = 2
    nx = 1
    ny = 2
  []
  [block1]
    type = SubdomainBoundingBoxGenerator
    input = 'msh'
    bottom_left = '0 0 0'
    top_right = '1 1 0'
    block_id = 1
    block_name = 'block1'
  []
  [block2]
    type = SubdomainBoundingBoxGenerator
    input = 'block1'
    bottom_left = '0 1 0'
    top_right = '1 2 0'
    block_id = 2
    block_name = 'block2'
  []
  [split]
    type = BreakMeshByBlockGenerator
    input = block2
  []
  [top_node]
    type = ExtraNodesetGenerator
    coord = '0 2 0'
    input = split
    new_boundary = top_node
  []
  [bottom_node]
    type = ExtraNodesetGenerator
    coord = '0 0 0'
    input = top_node
    new_boundary = bottom_node
  []
[]

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Modules]
  [TensorMechanics]
    [Master]
      generate_output = 'stress_yy'
      [all]
        strain = FINITE
        add_variables = true
        use_automatic_differentiation = true
        decomposition_method = TaylorExpansion
      []
    []
  []
[]

[BCs]
  [fix_x]
    type = DirichletBC
    preset = true
    value = 0.0
    boundary = bottom_node
    variable = disp_x
  []

  [fix_top]
    type = DirichletBC
    preset = true
    boundary = top
    variable = disp_x
    value = 0
  []

  [top]
    type = FunctionDirichletBC
    boundary = top
    variable = disp_y
    function = 'if(t<=0.3,t,if(t<=0.6,0.3-(t-0.3),0.6-t))'
    preset = true
  []

  [bottom]
    type = DirichletBC
    boundary = bottom
    variable = disp_y
    value = 0
    preset = true
  []
[]

[Modules/TensorMechanics/CohesiveZoneMaster]
  [czm_ik]
    boundary = 'interface'
  []
[]

[Materials]
  [stress]
    type = ADComputeFiniteStrainElasticStress
  []
  [elasticity_tensor]
    type = ADComputeElasticityTensor
    fill_method = symmetric9
    C_ijkl = '1.684e5 0.176e5 0.176e5 1.684e5 0.176e5 1.684e5 0.754e5 0.754e5 0.754e5'
  []
  [czm]
    type = BiLinearMixedModeTraction
    boundary = 'interface'
    penalty_stiffness = 1e6
    GI_c = 1e3
    GII_c = 1e2
    normal_strength = 1e4
    shear_strength = 1e3
    displacements = 'disp_x disp_y'
    eta = 2.2
    viscosity = 1e-3
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient

  solve_type = 'NEWTON'
  line_search = none

  petsc_options_iname = '-pc_type'
  petsc_options_value = 'lu'

  automatic_scaling = true

  l_max_its = 2
  l_tol = 1e-14
  nl_max_its = 30
  nl_rel_tol = 1e-50
  nl_abs_tol = 1e-15
  start_time = 0.0
  dt = 0.1
  end_time = 1.0
  dtmin = 0.1
[]

[Outputs]
  exodus = true
[]