Modeling of environmentally assisted intergranular crack propagation in polycrystals

Polycrystalline nickel‐based superalloys experience accelerated intergranular crack growth when exposed to dwell times in oxygen‐rich environments and a combination of high temperature and tensile mechanical loading. Increasing crack growth rates are observed for increasing amounts of environmental oxygen in a certain oxygen concentration range, while below and above that range crack growth rates remain approximately constant. A fully coupled chemo‐mechanical modeling framework accounting for the degradation of grain boundaries by oxygen has been presented by the authors. In this work, we expand the framework by a moving boundary condition to capture a realistic oxygen flux in grain boundary cracks for both edge cracks connected to the environment and interior cracks. In numerical simulation results, the behavior of the moving boundary condition is shown for intergranular crack propagation through a polycrystal subjected to cyclic loading. Finally, the capabilities of the modeling framework to qualitatively predict the dependence of the average crack growth rate on the environmental oxygen content, load level, and dwell time are evaluated and it is shown that predictions qualitatively agree with experimental observations for intergranular fracture.


INTRODUCTION
Nickel-based superalloys are often used in applications that require high-strength materials operating at high temperatures, such as turbine disks in jet engines.Under a combination of high temperatures, tensile mechanical loading and exposure to oxygen during significant dwell times, polycrystalline nickel-based superalloys are known to experience brittle intergranular failure.In these cases, the crack growth rate significantly increases in the presence of oxygen, possibly leading to accelerated failure.5][6] For both cases, this transition in crack growth rate coincides with a transition from transgranular crack growth (low oxygen levels/short or no dwell time) to intergranular crack growth (high oxygen levels/long dwell times); there is also a transition region with mixed trans-and intergranular crack growth when the crack growth rate starts to increase. 2,7Additionally it has been observed that the crack growth rate increases for increasing stress intensity factors for intergranular crack growth. 2,5,7 numerical modeling framework for describing the interaction of oxygen diffusion and mechanical material degradation has been suggested in Auth et al. 8 The framework includes a fully chemo-mechanically coupled cohesive law, describing the acceleration of mechanical degradation by the presence of oxygen as well as the acceleration of oxygen transport by the mechanical stress state.However, it does not yet provide a complete set-up for simulating crack propagation through polycrystals in a physically realistic manner.In particular, as cracks open up the outer boundary of the structure grows.Since oxygen is supplied form the environment, the domain where the boundary conditions on the oxygen field should be applied changes.Furthermore, the presented framework does not cover the computation of crack growth rates, where a particular problem is to find a criterion for determining the crack tip location.
In this work, we address these problems with the goal to evaluate crack propagation rates in cyclic loading.Especially, the chemo-mechanical modeling framework will be evaluated in its capability to correctly predict qualitative dependencies of crack growth rates on environmental conditions.For the scope of this manuscript, we focus on the dependence of crack growth rates on the environmental oxygen concentration, as experimentally investigated by for example, Molins et al. 2 Some approaches to boundary conditions on the chemical field in chemo-mechanically coupled problems can be found in the literature.In the context of hydrogen embrittlement, Serebrinsky et al. 9 have suggested a Dirichlet boundary condition on the crack flanks that is dependent on the hydrostatic stress.A similar approach has later been adopted by del Busto et al. 10 for cohesive zone modeling of hydrogen embrittlement.However, oxygen embrittlement differs from hydrogen embrittlement in that oxygen only diffuses along grain boundaries and not inside grains, 11 which makes this strategy less applicable to the present problem.In the context of oxygen assisted crack growth, Zhao 12 and Karabela et al. 13 have prescribed a time-dependent decrease in oxygen flux on the outer boundary of the structure, motivated by oxidation kinematics.This procedure, however, does not cover realistic oxygen concentrations in evolving cracks.
More recently, Kristensen et al. 14 have suggested penalizing the difference between the local and the environmental oxygen concentration in cracks in order to model hydrogen flow into edge cracks.This moving boundary condition has also been used by Golahmar et al. 15 for hydrogen assisted fatigue crack propagation.While the approach seems suitable for modeling the oxygen supply in edge cracks, its major drawback is that it introduces an oxygen source term and thus leads to unphysical behavior in the case of crack initiation inside the structure.
In this work, we develop a moving boundary condition for oxygen flow into edge cracks.The moving boundary condition is designed to handle the initiation of cracks on the domain boundaries as well as inside the structure.In order to detect cracked regions, a point-wise applicable, irreversible fracture criterion for cohesive zones is presented.Finally, crack propagation rates in experiments are often measured for quasi-static crack propagation controlled by the stress intensity factor.Therefore, an algorithm for prescribing the stress intensity factor in terms of force-controlled loading is incorporated into the modeling framework.

COUPLED CHEMO-MECHANICAL COHESIVE ZONE MODEL
In this section, a short overview over the employed coupled chemo-mechanical cohesive zone model is given.For a more detailed description, the reader is referred to Auth et al. 8

Kinematics
The displacement jump across the interfaces is given by the difference between the displacements at opposite sides of the interface u + and u − We perform a split of the displacement jump such that  = [ t , Δ n ].The displacement jumps tangential to the mid-plane of the interface  t and normal to it, Δ n , are where ñgb is the unit normal to the interface (compare Figure 1).Similarly, the gradient operator on the interfaces is defined as The oxygen flux in the grain boundaries j gb is given by where D is is the diffusion coefficient, c gb is the oxygen concentration in the mid-plane of the cohesive interface, R is the universal gas constant, T is the absolute temperature and  is the chemical potential for oxygen diffusion.

Balance equations
The balance equations in the bulk domain Ω and grain boundaries Γ s are the balance of linear momentum and conservation of oxide mass.The balance of linear momentum in Ω, for quasi-static conditions without body forces, is given as while on the grain boundaries Γ s it is given as where T + and T − are the traction vectors on opposite sides of the grain boundaries.The conservation of oxide mass in the grain boundaries Γ s , in the absence of a source term, gives

Chemo-mechanically coupled free energy
The free energy per unit area for the chemo-mechanically coupled interface material model Ψ Γ s is given by Auth et al. 8 as Therein,  0 is the reference chemical potential, N gb is the number of lattice sites per grain boundary area,  L = c gb ∕N gb is the lattice site occupancy, V O 2 is the partial molar volume of oxygen in metal, d O 2 is the environmental damage variable, K n and K t are cohesive stiffnesses in the normal and tangential directions and c gb 0 is a reference concentration.The stiffnesses K t and K n are chosen according to a modified version of the Xu-Needleman traction-separation law 16 with the irreversible damage unloading framework suggested by Kolluri et al. 17 The expressions for the stiffnesses are given by where Φ t and Φ n represent the tangential and normal work of separation and  t and  n are the characteristic tangential and normal separation jump.The damage variables are computed as Therein, the history variables for the maximum experienced jump  max and the maximum experienced oxygen concentration c gb max evolve as follows Furthermore, d O 2 ,max is the maximum environmental damage and c gb char is the characteristic oxygen concentration around which the environmental damage evolves.While this formulation allows for irreversible cyclic loading, additional modeling would be needed in order to incorporate fatigue characteristics such as the fatigue threshold or the endurance limit, compare for example, Reference 15.
The cohesive tractions T = [T t , T n ] and the oxygen flux in the grain boundaries j gb follow from the thermodynamical framework presented in Reference 8 together with the free energy in Equation ( 8) and the expression for the flux in Equation (4) as K n Δ n is the environmentally degraded normal traction.

MODELING OF CRACK PROPAGATION
The model presented in Auth et al. 8 describes the interaction of oxygen diffusion and the mechanical stress state ahead of the crack tip.The oxygen boundary conditions are prescribed on the outer boundary of the structure.In order to propagate cracks through the structure, it is important to discuss how oxygen boundary conditions should be handled in evolving edge cracks.Physically, the oxygen supply comes from the environment surrounding the structure.For cracks which open towards the domain boundary, we assume that oxygen flows in instantaneously compared to the time scale of oxygen flow in unbroken regions.Thus, the oxygen concentration in edge cracks should be equal to the environmental oxygen concentration up to the crack tip.Ahead of the crack tip, the oxygen transport is determined by the modeled multi-physics interaction and should not be influenced by the (moving) boundary conditions.Cracks are in general likely to initiate at environmentally exposed domain boundaries.However, on the microscale it is important to consider that (micro-)cracks can initiate inside the structure as well, for example, due to material defects.These cracks might have significant impact on for example, the crack growth pattern and the crack growth rate.To elaborate, consider two different fracture scenarios, as shown in Figure 2: in the first scenario (green/left) the structure breaks at a boundary, in the second scenario (orange/right) the structure breaks in the interior (e.g., due to a material weakness).From a physical perspective, the green crack should fill up with oxygen from the environment, contrarily the orange crack should not experience the environmental oxygen concentration.However, a significantly increased diffusivity can be assumed in the orange region compared to the undamaged grain boundaries.

Fracture criterion
In order to determine which parts of the grain boundaries are broken, a criterion for determining the fracture state of each material point is needed.This criterion will then be used for two purposes: (1) determine the domain where the moving boundary condition should be active, that is the fractured grain boundary domain Γ f .(2) Compute the crack length in order to determine the crack propagation rate and to control the loading via the stress intensity factor.In this section, a scalar, point-wise and irreversible fracture criterion is developed for the coupled cohesive law.The concept could in general be applied to other cohesive laws, but it is here presented for the specific formulation from Reference 8.
A fully broken cohesive zone does not produce any traction response, except for the case of contact between the interfaces (i.e., for compression).It is therefore natural to consider a traction based fracture criterion.We define a grain boundary point as broken once the remaining normal strength T r falls below a limit value T f .The remaining normal strength is the maximum possible traction response that the cohesive law can yield in the normal direction, considering the current damage state.
Consider the normal cohesive law, as shown in Figure 3.As long as Δ n <  n , the remaining normal strength T r is clearly the maximum of the traction-separation curve, that is  max .Once the normal jump Δ n has exceeded the characteristic normal jump  n , the remaining strength corresponds to the normal traction at the maximum experienced normal jump Δ n,max .The remaining normal strength T r can mathematically be formulated based on the environmentally degraded normal traction Tn = Tn For a given limit traction value for fracture T f , the broken regions Γ f of the grain boundaries Γ s are defined as Remark: The presented cohesive law has a strong coupling between the normal and tangential directions.In particular, a tangential separation jump leads to a decrease of both the normal and the tangential traction.For this reason, we decide to solely use the remaining normal strength in order to determine the interfaces fracture state.In general, the concept of remaining strength can similarly be applied to the tangential plane.In that case, the remaining normal and tangential strength both need to be accounted for in the scalar fracture criterion.

Penalty approach for the moving boundary condition
Different approaches can be imagined for applying the environmental oxygen concentration to the regions identified as broken, for example the modification of the Dirichlet boundary condition based on an active set search, a Lagrange multiplier approach or a penalty method.For the present model and its application, there is typically a large difference in oxygen scales: the environmental oxygen concentration is often much higher than the characteristic oxygen concentration needed for damaging the grain boundaries.As the oxygen in the crack represents the supply of oxygen at the crack tip, it is crucial for correct model behavior to capture this large concentration difference, although it is not crucial to capture the exact oxygen concentration in the crack.
Taking the previous points into consideration, we choose a penalty approach for ensuring that the environmental oxygen concentration is obtained in the edge cracks.As penalized quantity, the spatial gradient of concentration in the broken regions Γ f is selected, such that the penalty terms aims at fulfilling In combination with a Dirichlet boundary condition on the outer boundary of the domain Γ s , this approach enforces the environmental oxygen concentrations in edge cracks by suppressing spatial changes in the oxygen concentration.
For cracks inside the domain, it however only enforces a homogeneous oxygen level in the crack, which depends on the oxygen distribution in the surrounding grain boundaries.We implement this by adding a penalty term W bc to the backward Euler integrated weak form of the mass balance Equation (7) (see Reference 8 and Appendix B for details): where c mid is the test function for the concentration field.Taking the fracture criterion formulated in Section 3.1 into account, a penalty term for fulfilling Equation ( 16) can be expressed as where p bc is the penalty coefficient.Notice that this is equivalent to an additional diffusive flux j gb bc = D bc  gb c gb , with a non-constant diffusion coefficient D bc = p bc (T f − T r ).The contribution of the moving boundary condition to the weak form of the mass balance then reads Appendix B shows the derivation of this weak form in the context of the coupled cohesive element formulation presented in Reference 8.However, from a numerical perspective, Equation ( 18) is prone to cause problems.The Heaviside function is not continuously differentiable, which often poses problems for Newton-type solvers.Therefore, it is replaced by a regularized Heaviside function, that is continuously differentiable.The following regularized Heaviside function has been employed for the purpose of this work but it is generally possible to use other sigmoid functions.Additionally, in order to instantaneously fill edge cracks with oxygen the oxygen flux resulting from the penalty term has to be much larger than the other oxygen fluxes in the model.As a consequence, the moving boundary condition results in a significant jump of the oxygen flux at crack tips, causing steep jumps in the concentration field that are often unresolvable by mesh sizes sufficient to capture the chemo-mechanical fracture behavior ahead of the crack tip.The possibility to reduce this jump by the regularized Heaviside function is limited, as the remaining strength T r in the present model drastically increases in a narrow region ahead of the crack tip.It is further undesirable to introduce increased mesh requirements by the moving boundary condition.
For this reason, a smooth decrease in flux caused by the moving boundary condition over a wider region is sought for when approaching the crack tip.Opposed to the remaining strength, the displacement jump typically increases smoothly behind the crack tip and is therefore much better suited for smoothening the oxygen flow caused by the moving boundary condition.An additional advantage of using the normal displacement jump for regularization is that the additional diffusivity is decreased for small jumps and vanishes below a certain level, that is, no increased diffusivity occurs when a broken grain boundary is compressed.
Hence, the diffusion coefficient D bc is proposed to have the following form wherein T reg f and  reg f are the regularization parameters for the respective Heaviside functions.In the smoothening term,  f sets a level of the normal jump around which the regularization occurs.The regularization level  f should be chosen such that it is representative of the displacement jump at the crack tip.
The coupled cohesive zone model represents the material behavior ahead of the crack tip.Contrarily, the moving boundary condition describes oxygen flow behind the crack tip, that is, where the material is broken.In particular, the moving boundary condition should not impact the oxygen transport mechanisms ahead of the crack tip and instead be a strictly trailing phenomenon to material failure.This property was originally guaranteed by the fracture criterion, but has been lost in the process of regularization.Practically this means that ahead of the crack tip, the diffusion coefficient related to the moving boundary condition D bc should be small compared to the diffusion coefficient related to Fickian diffusion and traction-assisted oxygen flux D.
For discussing how to fulfill this requirement, we only take the first part of Equation ( 21 Equation ( 15) the crack tip is located at T r = T f and material points ahead of the crack tip are defined by T r > T f .We now introduce two additional parameters: a limit remaining strength T * r and an associated tolerated diffusivity D * tol .For all points where the remaining strength T r is larger than the limit remaining strength T * r , it is required that the diffusivity D bc is smaller than the tolerated diffusivity This behavior is achieved by introducing a modified regularized Heaviside function  mod reg such that where g(x) is a function that fulfills the following conditions: A piece-wise quadratic Ansatz is chosen for g(x) where it results from conditions 1 and 2 that b = 0 and c = 1.The remaining factor a can be determined iteratively based on condition 3 for given model parameters.Employing the modified Heaviside function, the diffusivity of the moving boundary condition finally becomes Note that the computational implementation here is much facilitated by the use of automatic differentiation.The solid line in Figure 4 shows the suggested modified Heaviside function.It is identical to the original regularized Heaviside function for remaining strength below the fracture traction T f and decreases to D * tol ∕p bc for the limit remaining strength T * r .Thus, by accepting a small transitioning regime between T f and T * r it is possible to choose large values for the penalty factor p bc without impacting the separation between pre-and post-fracture models while keeping a smooth numerical behavior.

NUMERICAL EXPERIMENTS
Numerical experiments are conducted to show the behavior of the proposed moving boundary condition and to evaluate crack growth rates in a polycrystal.There are two major aspects that are investigated here: (1) evaluation of the moving boundary condition regarding the requirements described in Section 3. ( 2) Testing if the chemo-mechanical model presented in Reference 8 combined with the presented moving boundary condition can produce crack growth rate patterns similar to those observed experimentally. 2or all simulations, a model problem with 15 grains with an average grain size of 150 m is used.The constitutive behavior of the grains is modeled by the crystal (visco)plasticity model and parameter values presented in Reference 8.The crystallographic directions are given in Appendix A. The mesh of the model problem and its boundary conditions are represented in Figure 5.The region marked by the green solid line is pre-fractured by prescribing an initial maximum normal jump and an initial oxygen concentration.
The polycrystal is loaded in traction controlled cyclic loading.The load cycles are defined by the stress-intensity factor and the prescribed traction is computed as described in Section 4.1.The stress-intensity factor is cycled in trapezoidal cycles with the load ratio R = 0.1.Time steps of 0.5 s during loading/unloading and of 2 s during the dwell time are used.The applied model parameters, as well as details on the initial conditions are given in Appendix A.
The finite element code is written in Julia, 18 relying on the tensor calculus toolbox Tensors.jl 19and the finite element toolbox Ferrite.jl. 20Grain structures are generated by Neper 21 and meshing is done with Gmsh. 22After meshing, Neper is again used for inserting cohesive elements along the grain boundaries.For result visualization, Makie.jl 23has been used.

Load control algorithm
The present numerical framework requires stable crack propagation.For achieving this we choose to prescribe an approximate stress intensity factor K I .This choice is in particular motivated by the analogy to experiments conducted by for example, Molins et al., 2 Moverare and Gustafsson, 4 Viskari et al. 24 The load is controlled in an explicit manner, that is, the load for the current time step n + 1 is computed based on the crack length of the previous time step n, that is, n a.
For a crack length a in a structure of width b the prescribed traction t * is based on linear fracture mechanics under the assumption of an edge crack in a plate subjected to uni-axial stress ) . (26) The crack length perpendicular to the loading direction is computed as follows: 1. Determine the fracture state of all material points according to Equation (15).
2. Starting at the pre-crack, find all crack paths in the structure.
3. Compute the crack lengths of all crack paths projected to the direction perpendicular to the loading direction.The longest crack path is chosen as the main fracture path and its projected length is used as the crack length a.
It should be pointed out that the control of the stress intensity factor K I is an approximation here, as several requirements for the above given analytical formula do not apply.In particular the material is not linear elastic and the crack growth direction is only approximately perpendicular to the loading direction (showing mixed mode crack growth).However, the procedure proves to work well for maintaining stable crack growth in practice.

Moving boundary condition
We will first examine the model problem, compare Figure 5, exposed to an environmental oxygen concentration of c env = c gb env ∕h = 2 nmol∕mm 3 , where h is the cohesive zone thickness, and a maximum stress intensity factor K I of 480 MPa √ mm.The load cycles consist of 10 s loading-300 s dwell time-10 s unloading.Figure 6A shows three snapshots of the von Mises stresses  vM in the structure during crack propagation.A pronounced stress concentration is observable around the crack tip, in the later time steps a second stress concentration develops on the upper edge of the structure due to bending.In Figure 6B, the oxygen concentrations along the projected crack coordinate s proj are displayed for the same time steps as the von Mises stresses are shown in A. The projected crack coordinate is a projection of the local crack directions to the direction normal to the loading direction, in this case x and s proj coincide.As the crack advances, the broken grain boundaries experience the environmental oxygen concentration.When approaching the crack tip, the oxygen concentration smoothly decreases.
The three different oxygen fluxes that transport oxygen into the structure are shown in Figure 6C.Oxygen can be transported by Fickian diffusion and by normal traction gradients ahead of the crack tip; the respective fluxes are here called j gb chem and j gb mech .Figure 6C shows the oxygen fluxes in the material point marked by a red cross in Figure 6A  the material point is broken and does not interact with the transport of oxygen while the material is degrading.It can further be observed that the traction-assisted oxygen flux j gb mech is dominant over the Fickian flux j gb chem during the initial oxygen transport.After material fracture, oxygen transport due to the moving boundary condition dominates the oxygen transport through the broken grain boundary.

Interior crack initiation
The presented moving boundary condition is further able to handle cracks that initiate inside the structure and are initially unconnected to the oxygen supply.In order to demonstrate this capability, a material weakness in the region marked by the orange line in Figure 5 is introduced.This is implemented by lowering the fracture energies in normal and tangential direction to 0.5 % of their original value, such that the region will break during the initial loading phase.The results are presented in Figure 7, the layout and the chosen material point are the same as in Figure 6.Notice that the weakened region surrounds the material point which is marked by the red cross.In Figure 7A, the crack in the interior of the structure can be recognized by the lower von Mises stresses above and below it, as well as the stress concentrations left and right to it (compare the first two time points).There are two noteworthy results in Figure 7B,C: (1) neither oxygen nor oxygen flux are present in the material point until the crack tip reaches it.This is opposite to what would happen in a penalty approach that penalizes c gb env − c gb , as for example, presented in Reference 14, where this scenario would lead to an oxygen source in the crack.(2) Once oxygen reaches the broken region in the interior, the oxygen flux caused by the moving boundary condition immediately increases.No oxygen flux due to mechanical interaction occurs, as the material point is already broken and cannot sustain any load.This means that oxygen is quickly transported through cracked regions once it reaches their boundary.This behavior is in accordance with the expected physical behavior upon crack initiation in the interior of the structure.In accordance with experimental results, crack growth rates are increasing faster at lower oxygen levels than at higher oxygen levels.

Crack growth rate during dwell time
avoid significant influence of crack branching that happens in some load cases at the grain boundary intersection located at s proj ≈ 0.25 mm.Analogously to experimental results presented in Reference 2, Figure 8 shows the obtained average crack growth rates per load cycle against the environmental oxygen concentrations in a double logarithmic diagram for different load levels.As the simulations are conducted on a microscale while the experiments are conducted on the macroscale, results are compared on the basis of capturing similar effects only.It is further assumed that the environmental influence on the crack growth rate is much more pronounced in this range of oxygen concentrations than fatigue effects.In accordance with the reported experimental results, we observe a larger increase of crack growth rate at lower oxygen concentrations than at higher oxygen concentrations.Higher load levels result in faster crack growth.Shorter dwell times lead to significantly less crack propagation per loading cycle compared to longer dwell times.

CONCLUDING REMARKS
From the outset of a fully coupled chemo-mechanical cohesive zone model, presented in Reference 8, we have in this article proposed a moving boundary condition for the chemical field in the grain boundaries.This formulation ensures physically meaningful oxygen concentrations in cracks, while allowing for crack initiation inside and on the boundary of the polycrystal.In conjunction with the moving boundary condition, a fracture criterion for the cohesive law has been suggested to identify the crack tip location.The fracture criterion further allows the use of an algorithm for controlling the load level for a given stress intensity factor while accounting for the crack length.Numerical experiments have shown that the model suggested in Reference 8, in combination with the proposed moving boundary condition for controlled stress intensity factors, can be used to simulate crack propagation during cyclic loading through polycrystals.The environmental oxygen concentration is reached in edge cracks and follows the crack tip upon propagation.It has been demonstrated that crack initiation inside the grain structure is handled physically correct by the moving boundary condition.The diffusivity in such an interior crack is increased; however, there is no direct impact of the moving boundary condition until the oxygen field reaches the crack.The numerical experiments have further shown that the original model behavior ahead of the crack tip is not disturbed by adding the moving boundary condition.
Finally, we have demonstrated by a series of cyclic loading simulations on a polycrystal that the crack rate pattern observed in intergranular fracture can be reproduced in a qualitatively correct way.In particular, the model qualitatively depicts the increase in average crack growth rate with increasing load level, dwell time, and oxygen levels.Faster increase of the average crack growth rate is predicted for lower than for higher environmental oxygen levels, thus capturing the saturation behavior seen in experiments.The range of constant crack growth rates at low oxygen concentrations, that is observed in experiments, is outside of the application range of this intergranular fracture model.

F I G U R E A1
Numbering of the grains, compare the crystal orientations in Table A2.

APPENDIX B. COUPLED COHESIVE ELEMENT FORMULATION
Similarly to the derivations in Appendix B of Reference 8, the weak form of the moving boundary condition can be split into two uncoupled equations by introducing a (small) cohesive zone thickness h assuming linear variation of the concentration perpendicular to the interfaces.In a local coordinate system  defined by the in-plane axes  1 and  2 and an axis  3 that is perpendicular to the interfaces, the concentration in 3D space can then be expressed as The in-plane concentration in the interfaces c gb is defined as Additionally to the gradient operator on the interface  gb , a gradient operator perpendicular to the interface  ⊥ is introduced such that  =  gb +  ⊥ .Starting from a general 3D version of the weak form of the moving boundary condition we insert the additive splits of the gradient operator  and the concentration c (compare Equation (B1)), as well as an analogous split of the variation of the concentration c.By integrating over the  3 -coordinate, we obtain

F I G U R E 2
Different fracture scenarios: The green edge crack (connected to the domain boundary) should experience the environmental oxygen concentration.The interior crack (inside the structure) should not be filled with oxygen.

F I G U R E 3
Loading/unloading curve under pure normal loading.Before reaching a normalized jump of Δ n ∕ n = 1, the remaining strength T r corresponds to the maximum strength of the cohesive law (compare point 1).For larger normal jumps, the remaining normal strength corresponds to the current normal traction upon loading (compare point 2) and to the traction at the beginning of unloading upon unloading/reloading (compare point 3).
) into account.The dashed line in Figure 4 exemplary shows  reg ( T f − T r , T reg r ) for increasing remaining strength T r .Remember that according to F I G U R E 4 Original and modified regularized Heaviside functions  reg and  mod reg , respectively.Function values of the modified Heaviside function are guaranteed to be smaller than D * tol ∕ p bc for remaining strength above the limit remaining strength T * r .This ensures that the moving boundary condition does not interfere with the model behavior ahead of the crack tip.

E 5
Model problem for numerical experiments.The structure is loaded under force control with a prescribed traction t * .An environmental oxygen concentration c env is prescribed on the left boundary of the marked grain boundary.The structure is pre-cracked by prescribing an initial maximum normal jump and an initial oxygen concentration along the solid green line.The region marked by the solid orange line inside the structure is weakened in some simulations in order to provoke a crack inside the structure.

F
I G U R E 6 (A) Evolution of the crack during the dwell time, displacements are magnified by factor 2. The red cross marks the material point represented in (C).The results are shown for the end of the hold period in the respective loading cycles.(B) Oxygen distribution along the crack at different time points.The projected coordinate along the crack s proj corresponds to the x-coordinate in this case.(C) Oxygen fluxes in the material point marked in (A).Notice that the fluxes are represented on a semi-logarithmic axis in order to display their size in relation to each other better.The time points used in (A) and (B) are marked by vertical lines.
over time.Notice the semi-logarithmic axes that are used in order to display all three fluxes despite their difference in magnitude.It can be seen that the chemical and mechanical oxygen fluxes j gb chem and j gb mech occur first, before the flux caused by the moving boundary condition j gb bc starts to increase.The flux caused by the moving boundary condition j gb bc occurs once F I G U R E 7 (A) Evolution of the crack during the dwell time, displacements are magnified by factor 2. The structure has been precracked around the red cross, as reflected by the lower von Mises stresses in that region.The red cross marks the material point represented in (C).The results are shown for the end of the hold period in the respective loading cycles.(B) Oxygen distribution along the crack at different time points.The projected coordinate along the crack s proj corresponds to the x-coordinate in this case.(C) Oxygen fluxes in the material point marked in (A).Notice that the fluxes are represented on a semi-logarithmic axis in order to display their size in relation to each other better.The time points used in (A) and (B) are marked by vertical lines.

Finally, we examine
crack growth rates through the model problem (without interior material weakness) for varying environmental oxygen levels.Cracks advance differently fast along the different grain boundaries, thus the average crack growth rate is determined based on a fixed spatial region.Here, projected coordinates s proj ∈ [0.08 mm, 0.197 mm] are used, in order to allow the crack to propagate away from the regions heavily influenced by the initial conditions and to F I G U R E 8 Crack growth rates for different maximum stress intensity factors K and increasing environmental oxygen concentrations.
c() = c mid ( 1 ,  2 ) + Δc( 1 , 2 )  h .(B1)Therein, the the average concentration c mid and the concentration jump Δc are defined by the concentrations at opposite sides of the interface c + and c − such thatc mid = 1 2 (c + + c − ) and Δc = c + − c − .(B2) 8he displacement field u is defined on the bulk domain Ω and the grain boundary (surface) domain Γ s , while the concentration field c is defined only on the grain boundaries Γ s .The normals to the outer boundaries of the domain are denoted n for the bulk boundary Ω and n gb for the normal to the outer boundary to the surface domain Γ s .The normal to the grain boundary domain Γ s is denoted ñgb .Figure originally published in Auth et al.8 Material parameter values employed for the numerical experiments.The values have been chosen such that qualitatively reasonable model behavior is obtained.Parameter values for the crystal plasticity model remain unchanged compared to Auth et al. 8 TA B L E A2 Rodrigues vectors employed in the polycrystalline example, compare the grain numbers to Figure A1. Note: ∫ Γ s D bc (  gb ( h c mid ) ⋅  gb c mid + h 12  gb Δc ⋅  gb Δc + Since the variations c mid and Δc are arbitrary, this can be split into two uncoupled weak forms W bc,gb = ∫ Γ s D bc  gb c gb ⋅  gb c mid dA, (B6)