Coupled Phase field with Mechanic tensor #20484

Judy-22 · 2022-03-07T17:42:08Z

Judy-22
Mar 7, 2022

Hi Moose community,

I am trying to model the coupling between phase fields and mechanics, using the material properties depending on the phase variables, but it converges. Below is the error shown.

0 Nonlinear |R| = 1.026418e+04
      0 Linear |R| = 1.026418e+04
  Linear solve did not converge due to DIVERGED_PC_FAILED iterations 0
                 PC failed due to SUBPC_ERROR 
Nonlinear solve did not converge due to DIVERGED_LINE_SEARCH iterations 0
 Solve Did NOT Converge!
  Finished Solving                                                                       [ 10.92 s] [    0 MB]
Aborting as solve did not converge
Solve failed, cutting timestep.

The input file:

[GlobalParams]        # Automatically define order parameters 
  var_name_base = eta
  op_num = 2
  displacements = 'disp_x disp_y'
[]
[Mesh]                # Generate mesh (7/3D X 4/3D) 
  type = GeneratedMesh
  dim = 2                                
  nx = 70             # n°s element in x. 
  ny = 40             # n°s element in y.
  xmin = 0  
  xmax = 70.0      # µm
  ymin = 0
  ymax = 40.0      # µm
  elem_type = QUAD4
  #uniform_refine = 2
[]

[Variables]
  [./c]               
    order = FIRST
    family = LAGRANGE
  [../]
  [./w]               
    order = FIRST
    family = LAGRANGE
  [../]
  [./PolycrystalVariables] 
  [../]
[]

[Kernels]
  [./c_res]
    type = SplitCHParsed
    variable = c
    f_name = F
    kappa_name = kappa_c
    args = 'eta0 eta1'
    w = w
  [../]
  [./w_res]
    type = SplitCHWRes
    variable = w
    mob_name = M
  [../]
  [./time]
    type = CoupledTimeDerivative
    variable = w
    v = c
  [../]
  [./RigidBodyMultiKernel]
    c = c
    f_name = F
    mob_name = L
    kappa_name = kappa_eta
    op_num = 2
    var_name_base = eta
    grain_force = grain_force
    grain_volumes = grain_volumes
    grain_tracker_object = grain_center
  [../]
  [./CH_RBM]
    type = MultiGrainRigidBodyMotion
    variable = w
    c = c
    op_num = 2
    var_name_base = eta
    grain_force = grain_force
    grain_volumes = grain_volumes
    grain_tracker_object = grain_center
  [../]
  [./vadv_eta0]
    type = SingleGrainRigidBodyMotion
    variable = eta0
    c = c
    v = 'eta0 eta1'
    grain_force = grain_force
    grain_tracker_object = grain_center
    grain_volumes = grain_volumes
    op_index = 0
  [../]
  [./vadv_eta1]
    type = SingleGrainRigidBodyMotion
    variable = eta1
    c = c
    v = 'eta0 eta1'
    op_index = 1
    grain_force = grain_force
    grain_tracker_object = grain_center
    grain_volumes = grain_volumes
  [../]
[]

[Modules/TensorMechanics/Master]
  [./all]
    add_variables = true
  [../]
[]

[Materials]
  [./Mat_prop]
    type = GenericConstantMaterial
    prop_names  = 'Es    Egb    Iw     F⁰       l⁰       t⁰
                   Dv0   Ds0   Dgb0    Qvol      Qsuf    Qgb
                   KB    Vm     Mgb    Qm        wGB     T '
    
    prop_values = '2.41        1.06      3.0e-6     1.0e-06  1.0e-06   1.0
                   2.0e-04     0.4       0.127e-04  2.6      2.21      0.58
                   8.6173e-05  1.21e-29  5.53e-08   0.171    3.0e-06   600'
  [../]
  [./Kc]
    type = ParsedMaterial
    f_name  = kappa_c
    function = '(l⁰/F⁰)*(3/4)*Iw*(2*Es - Egb)'
    material_property_names = 'Iw Es Egb F⁰ l⁰'
  [../]
  [./Keta]
    type = ParsedMaterial
    f_name  = kappa_eta
    function = '(l⁰/F⁰)*(3/4)*Iw*Egb'
    material_property_names = 'Iw Egb F⁰ l⁰'
  [../]
  [./A]
    type = ParsedMaterial
    f_name  = A
    function = '(l⁰^3/F⁰)*(12*Es - 7*Egb)/Iw'
    material_property_names = 'Iw Es Egb F⁰ l⁰'
  [../]
  [./B]
    type = ParsedMaterial
    f_name  = B
    function = '(l⁰^3/F⁰)*Egb/Iw'
    material_property_names = 'Iw Egb F⁰ l⁰'
  [../]
  [./Dvol]
    type = ParsedMaterial
    f_name  = Dvol
    function = '(t⁰/l⁰^2)*Dv0*exp(-Qvol/(KB*T))'
    material_property_names = 'Dv0 Qvol KB T t⁰ F⁰ l⁰'
  [../]
  [./Dsuf]
    type = ParsedMaterial
    f_name  = Dsuf
    function = '(t⁰/l⁰^2)*Ds0*exp(-Qsuf/(KB*T))'
    material_property_names = 'Ds0 Qsuf KB T t⁰ F⁰ l⁰'
  [../]
  [./Dgb]
    type = ParsedMaterial
    f_name  = Dgb
    function = '(t⁰/l⁰^2)*Dgb0*exp(-Qgb/(KB*T))'
    material_property_names = 'Dgb0 Qgb KB T t⁰ F⁰ l⁰'
  [../]
  [./chemical_energy]            
    type = DerivativeParsedMaterial
    f_name = Fc
    args = 'c eta0 eta1'          # Must be changed as op_num changes.
    function = 'A*c^2*(1-c)^2 + B*(c^2+6*(1-c)*(eta0^2+eta1^2)-4*(2-c)*(eta0^3+eta1^3)
                + 3*(eta0^2+eta1^2)^2)'
    material_property_names = 'A  B'
    derivative_order = 2
  [../]
  [./CHmobility]   
    type = DerivativeParsedMaterial
    f_name = M
    args = 'c eta0 eta1'
    function = '(t⁰*F⁰/l⁰^3)*Vm/(1.3806e-23*T)*(Dvol*c^3*(10 - 15*c + 6*c^2) + Dsuf*(1 - c) + Dgb*(eta0*eta1))'
    material_property_names = 'Dvol Qvol Dvap  Dsuf Dgb  KB Vm T t⁰ F⁰ l⁰'
    derivative_order = 1
  [../]
  [./L]       
    type = ParsedMaterial
    f_name  = L
    function = '(t⁰*F⁰/l⁰^3)*(4/3)*(Mgb/wGB)*exp(-Qm/(KB*T))'
    material_property_names = 'Mgb Qm KB T wGB t⁰ F⁰ l⁰'

  [../]
  [./force_density]  
    type = ForceDensityMaterial
    c = c
    etas ='eta0 eta1'
  [../]
  
  [./YoungsModulus]
    type = ParsedMaterial
    f_name  = E
    function = '(3*c^2 - 2*c^3)*210e09'
    #material_property_names = 'l⁰'
    args = 'c'
  [../]
  [./poissons_ratio]
    type = ParsedMaterial
    f_name  = poissons_ratio
    function = '(3*c^2 - 2*c^3)*0.3'
    args = 'c'
  [../]
  [./ElasticityTensor]
    type = ComputeVariableIsotropicElasticityTensor
    youngs_modulus = E # N/m²
    poissons_ratio = poissons_ratio
    args = c
  [../]
  [./stress]
    type = ComputeLinearElasticStress
  [../]
  [./elastic_free_energy]
    type = ElasticEnergyMaterial
    f_name = Fe
    block = 0
    args = 'c'
    derivative_order = 2
  [../]
  [./free_energy]
    type = DerivativeSumMaterial
    block = 0
    f_name = F
    sum_materials = 'Fc Fe'
    args = 'c eta0 eta1'
    derivative_order = 2
  [../]
[]

[BCs]
  [./bcs]
    #zero flux BC
    type = NeumannBC
    value = 0
    variable = c
    boundary = '0 1 2 3'
  [../]
  [./bottom_x]
    type = DirichletBC
    variable = disp_x
    boundary = bottom
    value = 0
  [../]
  [./bottom_y]
    type = DirichletBC
    variable = disp_y
    boundary = bottom
    value = 0
  [../]
  [./Pressure]
    [./right]
      boundary = 'right left top'
      function = 10e+3*t
    [../]
  [../]
[]

[AuxVariables]
  [./bnds]
  [../]
  [./free_energy]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./df00]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./df01]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./df10]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./df11]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./unique_grains]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./var_indices]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./centroids]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./vadvx]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./vadvy]
    order = CONSTANT
    family = MONOMIAL
  [../]
[]

[AuxKernels]
  [./bnds]
    type = BndsCalcAux
    variable = bnds
    var_name_base = eta
    op_num = 2
    v = 'eta0 eta1'
  [../]
  [./energy_density]
    type = TotalFreeEnergy
    variable = free_energy
    f_name = F
    kappa_names = 'kappa_c kappa_eta kappa_eta'
    interfacial_vars = 'c eta0 eta1'
  [../]
  [./df01]
    type = MaterialStdVectorRealGradientAux
    variable = df01
    index = 0
    component = 1
    property = force_density
  [../]
  [./df11]
    type = MaterialStdVectorRealGradientAux
    variable = df11
    index = 1
    component = 1
    property = force_density
  [../]
  [./df00]
    type = MaterialStdVectorRealGradientAux
    variable = df00
    index = 0
    component = 0
    property = force_density
  [../]
  [./df10]
    type = MaterialStdVectorRealGradientAux
    variable = df10
    index = 1
    component = 0
    property = force_density
  [../]
  [./unique_grains]
    type = FeatureFloodCountAux
    variable = unique_grains
    flood_counter = grain_center
    field_display = UNIQUE_REGION
    execute_on = timestep_begin
  [../]
  [./var_indices]
    type = FeatureFloodCountAux
    variable = var_indices
    flood_counter = grain_center
    field_display = VARIABLE_COLORING
    execute_on = timestep_begin
  [../]
  [./centroids]
    type = FeatureFloodCountAux
    variable = centroids
    execute_on = timestep_begin
    field_display = CENTROID
    flood_counter = grain_center
  [../]
  [./vadv_x]
    type = GrainAdvectionAux
    component = x
    grain_tracker_object = grain_center
    grain_force = grain_force
    grain_volumes = grain_volumes
    variable = vadvx
  [../]
  [./vadv_y]
    type = GrainAdvectionAux
    component = y
    grain_tracker_object = grain_center
    grain_force = grain_force
    grain_volumes = grain_volumes
    variable = vadvy
  [../]
[]

[ICs]
  [./ic_eta0]
    type = SmoothCircleIC
    variable = eta0
    x1 = 20.0
    y1 = 20.0
    radius = 15.0
    invalue = 1.0
    outvalue = 0.0
    int_width = 3.0
  [../]
  [./IC_eta1]
    type = SmoothCircleIC
    variable = eta1
    x1 = 50.0
    y1 = 20.0
    radius = 15.0
    invalue = 1.0
    outvalue = 0.0
    int_width = 3.0
  [../]
  [./ic_c]
    type = SpecifiedSmoothCircleIC
    invalue = 1.0
    outvalue = 0.0
    int_width = 3.0
    x_positions = '20.0 50.0'
    y_positions = '20.0 20.0'
    z_positions = '0.0 0.0'
    radii       = '15.0 15.0'
    3D_spheres = false
    variable = c
    block = 0
  [../]
[]

[VectorPostprocessors]
  [./forces]
    type = GrainForcesPostprocessor
    grain_force = grain_force
  [../]
  [./grain_volumes]
    type = FeatureVolumeVectorPostprocessor
    flood_counter = grain_center
    execute_on = 'initial timestep_begin'
  [../]
[]

[Postprocessors]
  [./average_c]
    type = ElementAverageValue
    variable = c
    execute_on = 'initial timestep_end'
  [../]
[]

[UserObjects]
  [./grain_center]
    type = GrainTracker
    outputs = none
    compute_var_to_feature_map = true
    execute_on = 'initial timestep_begin'
  [../]
  [./grain_force]
    type = ComputeGrainForceAndTorque
    execute_on = 'linear nonlinear'
    grain_data = grain_center
    force_density = force_density
    c = c
    etas = 'eta0 eta1'
  [../]
[]

[Preconditioning]
  [./SMP]
    type = SMP
    full = true
  [../]
[]

[Executioner]
  type = Transient
  scheme = bdf2
  solve_type = NEWTON
  petsc_options_iname = '-pc_type -ksp_gmres_restart -sub_ksp_type -sub_pc_type -pc_asm_overlap'
  petsc_options_value = 'asm         31   preonly   lu      1'
  l_tol = 1e-05
  nl_max_its = 30
  l_max_its = 30
  nl_rel_tol = 1e-07
  nl_abs_tol = 1e-09
  start_time = 0.0
  end_time = 10
  dt = 0.5
  [./TimeStepper]
    type = IterationAdaptiveDT
    cutback_factor = 0.5
    dt = 0.5
    growth_factor = 2
    optimal_iterations = 8
  [../]
[]

[Outputs]
  exodus = true
[]

Answered by SudiptaBiswas

Mar 7, 2022

Judy, You need to non-dimensionalize your mechanical properties based on the scales you are using in the phase-field part. Otherwise, the chemical free energy and elastic energy can vary by several orders of magnitude. This could cause numerical issues and also provide incorrect evolution. You can also turn on automatic_scaling = true in the Executioner block if the residuals for different nonlinear variables are of different orders of magnitude.

View full answer

GiudGiud · 2022-03-07T19:42:08Z

GiudGiud
Mar 7, 2022
Collaborator

Tagging @SudiptaBiswas @dschwen for help

0 replies

SudiptaBiswas · 2022-03-07T19:58:56Z

SudiptaBiswas
Mar 7, 2022

Judy, You need to non-dimensionalize your mechanical properties based on the scales you are using in the phase-field part. Otherwise, the chemical free energy and elastic energy can vary by several orders of magnitude. This could cause numerical issues and also provide incorrect evolution. You can also turn on automatic_scaling = true in the Executioner block if the residuals for different nonlinear variables are of different orders of magnitude.

11 replies

lindsayad Mar 8, 2022
Maintainer

Fix in #20502

Judy-22 Mar 8, 2022
Author

@SudiptaBiswas

I did as you recommended, I removed the the rigid-body motion, and now it works fine for automatic_scale. But the solve do not converge.

Time Step 1, time = 0.01, dt = 0.01
Performing automatic scaling calculation

0 Nonlinear |R| = 2.113839e+00
|residual|_2 of individual variables:
c: 1.60534
w: 0
eta0: 0.97242
eta1: 0.97242
disp_x: 1.82826e-05
disp_y: 1.72026e-05
0 Linear |R| = 2.113839e+00
Linear solve did not converge due to DIVERGED_PC_FAILED iterations 0
PC failed due to SUBPC_ERROR

Judy-22 Mar 11, 2022
Author

@SudiptaBiswas

I did as you recommended, I removed the the rigid-body motion, and now it works fine for automatic_scale. But the solve do not converge.

Time Step 1, time = 0.01, dt = 0.01
Performing automatic scaling calculation

0 Nonlinear |R| = 2.113839e+00
|residual|_2 of individual variables:
c: 1.60534
w: 0
eta0: 0.97242
eta1: 0.97242
disp_x: 1.82826e-05
disp_y: 1.72026e-05
0 Linear |R| = 2.113839e+00
Linear solve did not converge due to DIVERGED_PC_FAILED iterations 0
PC failed due to SUBPC_ERROR

SudiptaBiswas Apr 4, 2022

@Judy-22 Did the automatic scaling fix help? Are your material properties for the mechanical and phase-field part of the model are of similar orders of magnitude? Could you please post the new input file so I can test it?
Also, here is some guidance about failed solves, https://mooseframework.inl.gov/application_usage/failed_solves.html

Judy-22 Apr 4, 2022
Author

[sintering_316L_wRB_mec.zip]
@SudiptaBiswas,
See attached the input file.

(https://github.com/idaholab/moose/files/8412147/sintering_316L_wRB_mec.zip)

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Coupled Phase field with Mechanic tensor #20484

{{title}}

{{editor}}'s edit

{{editor}}'s edit

Replies: 2 comments 11 replies

{{title}}

{{title}}

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

Select a reply

Coupled Phase field with Mechanic tensor #20484

Judy-22 Mar 7, 2022

Replies: 2 comments · 11 replies

GiudGiud Mar 7, 2022 Collaborator

SudiptaBiswas Mar 7, 2022

lindsayad Mar 8, 2022 Maintainer

Judy-22 Mar 8, 2022 Author

Judy-22 Mar 11, 2022 Author

SudiptaBiswas Apr 4, 2022

Judy-22 Apr 4, 2022 Author

Judy-22
Mar 7, 2022

Replies: 2 comments 11 replies

GiudGiud
Mar 7, 2022
Collaborator

SudiptaBiswas
Mar 7, 2022

lindsayad Mar 8, 2022
Maintainer

Judy-22 Mar 8, 2022
Author

Judy-22 Mar 11, 2022
Author

Judy-22 Apr 4, 2022
Author