Question on Phase Field Fracture simulation results of the first time step

[Lingfu-Liu]

I created a model coupled solid mechanics and phase field method to study 1D bar elongation with damage by gradually adding loads (specify strain). I notice a very strange phenomenon that the first time step can not creat any damage no matter how huge strain I specify.

My script is attached as below.

[Mesh]
type = GeneratedMesh
dim = 1
nx = 50
xmin = 0.0
xmax = 1.0
[]

[GlobalParams]
displacements = 'disp'
[]

[Variables]

[./disp] # displacement
order = FIRST
family = LAGRANGE
[../]

[./d] # Phase field
order = FIRST
family = LAGRANGE
[../]
[]

[Physics/SolidMechanics/QuasiStatic]
[./mech]
strain = SMALL
add_variables = true
generate_output = stress_xx
material_output_family = LAGRANGE
material_output_order = FIRST
[../]
[]

[Modules/PhaseField/Nonconserved]
[./d]
free_energy = F
kappa = kappa_op
mobility = L
[../]
[]

[Kernels]
[./couple]
type = PhaseFieldFractureMechanicsOffDiag
variable = disp
component = 0
c = d
[../]
[]

[BCs]
[./left_fixed]
type = DirichletBC
boundary = left
variable = disp
value = 0.0
[../]

[./xdisp]
type = FunctionDirichletBC
variable = disp
boundary = right
function = 't'
[../]
[]

[Materials]
[./elasticity_tensor]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 100e9 # Pa
poissons_ratio = 0
[../]

[./pfbulkmat]
type = GenericConstantMaterial
# gc_prop: Critical energy release rate (Energy required to create fractures per unit area)
# l: internal length scale (diffusion leng;pth paramter)
# visco: damping coefficient of phase field
# c0: paramter in crack geometrical function (a(d)=d^2 ==> c0 =2)
prop_names = 'gc_prop l visco c0'
prop_values = '1 0.1 1e-4 2'
[../]
[./define_mobility]
type = ParsedMaterial
property_name = L
material_property_names = 'visco'
expression = '1.0/visco'
[../]
[./define_kappa]
type = ParsedMaterial
property_name = kappa_op
material_property_names = 'gc_prop l c0'
expression = '2 * gc_prop * l / c0'
[../]
[./elastic]
type = ComputeLinearElasticPFFractureStress
c = d
E_name = 'elastic_energy'
D_name = 'degradation'
F_name = 'local_fracture_energy'
decomposition_type = strain_spectral
[../]

[./degradation]
type = DerivativeParsedMaterial
property_name = degradation
coupled_variables = 'd'
expression = (1-d)^p*(1-delta)+delta
constant_names = 'p delta'
constant_expressions = '2 1e-6'
derivative_order = 2
[../]
[./local_fracture_energy]
type = DerivativeParsedMaterial
property_name = local_fracture_energy
coupled_variables = d
material_property_names = 'gc_prop c0 l'
expression = 'gc_prop * d^2 / c0 / l'
derivative_order = 2
[../]
[./fracture_driving_energy]
type = DerivativeSumMaterial
property_name = F
coupled_variables = d
sum_materials = 'elastic_energy local_fracture_energy'
derivative_order = 2
[../]

[]

[Executioner]
type = Transient
dt = 1e-6
solve_type = PJFNK
end_time = 1e-5
[]

[Outputs]
exodus = true
[./csv]
type = CSV
[../]
[]

[Postprocessors]
[u(end)]
type = PointValue
point = '1 0 0'
variable = disp
execute_on = TIMESTEP_END
[]

[d_ave]
type = ElementAverageValue
variable = 'd'
execute_on = TIMESTEP_END
[]

[sigma_xx_ave]
type = ElementAverageValue
variable = 'stress_xx'
execute_on = TIMESTEP_END
[]
[]

Postprocessor Values:
(Table 1)
+----------------+----------------+----------------+----------------+
| time | d_ave | sigma_xx_ave | u(end) |
+----------------+----------------+----------------+----------------+
| 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 |
| 1.000000e-06 | 0.000000e+00 | 1.000000e+05 | 1.000000e-06 |
| 2.000000e-06 | 9.082643e-04 | 1.996369e+05 | 2.000000e-06 |
| 3.000000e-06 | 4.445888e-03 | 2.973384e+05 | 3.000000e-06 |
| 4.000000e-06 | 1.212433e-02 | 3.903593e+05 | 4.000000e-06 |
| 5.000000e-06 | 2.520100e-02 | 4.751166e+05 | 5.000000e-06 |
| 6.000000e-06 | 4.462309e-02 | 5.476471e+05 | 6.000000e-06 |
| 7.000000e-06 | 7.097100e-02 | 6.041665e+05 | 7.000000e-06 |
| 8.000000e-06 | 1.044134e-01 | 6.416605e+05 | 8.000000e-06 |
| 9.000000e-06 | 1.446850e-01 | 6.584077e+05 | 9.000000e-06 |
| 1.000000e-05 | 1.910964e-01 | 6.543253e+05 | 1.000000e-05 |
+----------------+----------------+----------------+----------------+

In above table, "time" is equavalent to "strain" specified to the system. I manually calculate the analitical solution of "strain" versus "d"(phase filed of damage) with the belowing implicite process.

I have the following results.

(Table 2)
strain | d_analitical
0.00E+00 | 0
1.00E-06 | 0.000908265
2.00E-06 | 0.004445892
3.00E-06 | 0.012124339
4.00E-06 | 0.025201021
5.00E-06 | 0.04462313
6.00E-06 | 0.070971065
7.00E-06 | 0.104413459
8.00E-06 | 0.144685103
9.00E-06 | 0.191096615
1.00E-05 | 0.242580513

By comparing results in Tables 1 and 2, you may find that if I delete the data of "d" at the second row and second column in Table 1, the rest data of "time" vs "d_ave" in Table 1 matches well with "strain" vs "d_analitical" in Table 2.

In summary, my question is why there is a zero damage (d=0) results at first time step. Even I specify a very huge strain (e.g., 1), "d" is always equal to 0 at first time step.

[GiudGiud]

Hello

The first result in the table (if the table comes from MOOSE console output) in not from the first time step, but rather the INITIAL step.
And by default, no object is computed on INITIAL.
To get computations on initial, you need to run the postprocessors / auxkernels / etc (anything involved in computing this value) on INITIAL.
this can be done by setting execute_on = 'INITIAL TIMESTEP_END' for example

Note that at t=0 (=strain = 0 here?), you expect to only see any damage?

EDIT: my bad. I misread the last table as the output from MOOSE

[Lingfu-Liu]

Thank you for your response.

I included execute_on = TIMESTEP_END in all postprocessors - referring to my scripts.

I need to clarify that the first time step means t=dt=1e-6=strain. I see the stress becomes none zero at end of the first time step, but the damage (d) keeps zero. Even I set the strain=dt=1, the value of d always keeps zero at the end of first time step. I don't know how it works.

I am thinking if the problem is about the solving processes. Does MOOSE apply a fully coupled solver, i.e. solve two variables (disp and d in my case) simultaneously with one Jacobin matrix? Or does MOOSE solve "disp" first while keeping "d" as constant in the [QuasiStatic] Block, then update "d" with the new "disp" in the [Nonconserved] Block, while keeping the "disp" as constant (input) and only solving "d"?

[GiudGiud]

Hello

I dont see anything in the input file that deviates from the fully coupled apprach. So it should solve all variables together

[Lingfu-Liu]

Hi,

Thank you for your answering.

I created a new test simulation by specifying initial displacement (u_i=1e-6*x) while keeping d_i=0. Then I add loads by specifing u_end = dt+1e-6. I found out the value of "d" also keeps at 0 at end of the first time step (as shown in the second row and second column of the following Postprocessor results), which means MOOSE doesn't run with my assumed non-fully coupled method.

Do you have any idea why phase field variable "d" is not updated in the first time step?

Postprocessor Values:
+----------------+----------------+----------------+----------------+
| time | d_ave | sigma_xx_ave | u(end) |
+----------------+----------------+----------------+----------------+
| 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 1.000000e-06 |
| 1.000000e-06 | 0.000000e+00 | 2.000000e+05 | 2.000000e-06 |
| 2.000000e-06 | 3.623185e-03 | 2.978300e+05 | 3.000000e-06 |
| 3.000000e-06 | 1.138249e-02 | 3.909458e+05 | 4.000000e-06 |
| 4.000000e-06 | 2.453626e-02 | 4.757648e+05 | 5.000000e-06 |
| 5.000000e-06 | 4.403221e-02 | 5.483247e+05 | 6.000000e-06 |
| 6.000000e-06 | 7.045086e-02 | 6.048432e+05 | 7.000000e-06 |
| 7.000000e-06 | 1.039607e-01 | 6.423093e+05 | 8.000000e-06 |
| 8.000000e-06 | 1.442961e-01 | 6.590065e+05 | 9.000000e-06 |
| 9.000000e-06 | 1.907671e-01 | 6.548582e+05 | 1.000000e-05 |
| 1.000000e-05 | 2.423059e-01 | 6.315109e+05 | 1.100000e-05 |
+----------------+----------------+----------------+----------------+

[GiudGiud]

I'm afraid not. @laagesen would you know why?

[laagesen]

sorry I'm afraid not, maybe @SudiptaBiswas or @jiangwen84 ?