AN OPEN-SOURCE FENICS IMPLEMENTATION OF A PHASE FIELD FRACTURE MODEL FOR MICROPOLAR CONTINUA

A micropolar phase field fracture model is implemented in an open source library FEniCS. This implementation is based on the theoretical study in Suh, H.S., Sun, W., and O’Connor, D. (under review) in which the resultant phase field model exhibits the consistent micropolar size effect in both elastic and damage regions identifiable via inverse problems for micropolar continua. By leveraging the automatic code generation technique in FEniCS, we provide a documentation of the source code expressed in a language very close to the mathematical expressions without comprising significant efficiency. This combination of generality and interpretability therefore enables us to provide a detailed walk-through that connects the implementation with the regularized damage theory for micropolar materials. By making the source code open source, the paper will provide an efficient development and educational tool for third-party verification and validation, as well as for future development of other higher order continuum damage models.

Since the size effect is stemmed from the microstructures and topological features, the classical phenomenological constitutive laws, which rely on the evolution of local internal variables to represent history, may not be sufficient (Fleck et al., 1994;Mota et al., 2013;Miehe et al., 2013;Lin et al., 2015). Furthermore, the size effect also makes the first-order homogenization not suitable, because the corresponding effective medium does not capture the highorder kinematics Na et al., 2019;Hu and Oskay, 2019). From the numerical perspective, rate-independent constitutive laws formulated at material points may also lead to spurious mesh-dependence when softening and/or damage leads to strain localization (Bazant et al., 1984;Song et al., 2008;Na and Sun, 2016). Consequently, regularization procedure is needed to ensure that the boundary value problem remains well-posed (Belytschko et al., 2013;Sun and Mota, 2014).
The phase field or other nonlocal fracture models [e.g., eigenfracture (Schmidt et al., 2009;Wang and Sun, 2017;Qinami et al., 2019) and thick level set (Moës et al., 2011;Cazes and Moës, 2015)] provides a simple treatment to introduce size effect for the damage mechanics. By employing the smooth implicit function to represent a sharp crack interface, this family of models regularizes the sharp crack surfaces by diffusive damage zones and therefore bypass the need to embed strong discontinuities (Bourdin et al., 2008;Kuhn and Müller, 2010;Miehe et al., 2010a,b;Hofacker and Miehe, 2013;Choo and Sun, 2018;Bryant and Sun, 2018). This implicit representation of a sharp interface also shows promising results in simulating crack branching and dynamic fractures (Borden et al., 2012). However, the introduction of the length scale to regularize the phase field also stimulates debates on the legitimacy and justification of the length scales and on the proper inverse problems for identifying the material parameters (Lorentz et al., 2012;Wu and Nguyen, 2018;Geelen et al., 2019;Wu et al., 2020). In a recent work [cf. Suh et al. (2020)], we introduce the phase field theory for cohesive fracture of micropolar continua. Our goal of this previous research was to provide a mean to introduce size effect directly from the higher order kinematics while introducing a mean to decouple the influence of the length scale parameter of the phase field.

Why Open Source?
While technically feasible, the implementation of a phase field model for micropolar continua can be a complex task. Even though the phase field fracture does not require embedding strong discontunities through remeshing or enriching basis function, it does introduce additional field variables that must be resolved. This coupled system of equations cannot be updated easily in a monolithic solver; hence, different designs of sequential solvers have been developed in the last several decades (Amor et al., 2009;Ambati et al., 2015). To extend a phase field model for micropolar continua, additional considerations must be given to solve the interpolated microrotation, displacement, and phase field properly.
In this work, our goal is to provide a detailed walk-through of a mixed finite element implementation of the micropolar phase field model with a code base called FEniCS that enables simulations with limited programming (Logg et al., 2012;Alnaes et al., 2015;Langtangen et al., 2016). By leveraging a domain-specific language embedded in Python, called unified form language, we specify the finite element discretization of the variational form of our models and design a sequential solver that provides incremental solutions in two-dimensional domains. By making this implementation completely accessible with an open-source license, we may potentially remove the technical barrier for practitioners and researchers who are not necessarily familiar with computer implementation but are interested at experimenting with the phase field fracture theory for micropolar continua. Consequently, this work will benefit users from the many built-in features of the code, including the many scientific computing tools embedded in FEniCS, while it enables them to adapt and extend the applications of the models to a broader class of problems. More importantly, by making the code transparent and accessible to the public, we provide an opportunity for third-party verification and validation such that an open standard can be established. By minimizing the waste of time to "reinvent the wheel," this open source implementation will also save time and resources. The metadata of the open-sourced program is summarized in Table 1.
As for notations and symbols, bold-faced and blackboard bold-faced letters denote tensors (including vectors, which are rank-one tensors); a centerdot denotes a single contraction of adjacent indices of two tensors (e.g., a · b = a i b i or c · d = c ij d jk ); a colon denotes a double contraction of adjacent indices of tensor of rank two or higher (e.g., C : ε = C ijkl ε kl ); the ⊗ denotes a juxtaposition of two vectors (e.g., a ⊗ b = a i b j ) or two symmetric second-order tensors [e.g., (α ⊗ β) ijkl = α ij β kl ]. We also define identity tensors: I = δ ij , I = δ ik δ jl andĪ = δ il δ jk where δ ij is the Kronecker delta. As for sign conventions, unless specified, the directions of the tensile stress and dilative pressure are considered as positive.

FENICS IMPLEMENTATION OF MICROPOLAR PHASE FIELD MODEL
This section outlines the necessary steps that implements the length-scale-insensitive phase field model for cohesive fracture in micropolar continua, based on our theoretical [presented in Suh et al. (2020)]. Our focus is on how to leverage the Unifed Form Language (UFL), a domain-specific language, to implement the variational form of a boundary value problem such that the information flow of the codes stay close to the mathematical notation as much as possible.
Throughout this section, we present a step-by-step mathematical derivation, where at each step we provide the corresponding part of FEniCS implementation. Since all the scripts we developed are available online, we only illustrate and explain the code blocks that have a direct one-to-one relationship with the derivation steps.
We start with reviewing the kinematic and constitutive relations for size-dependent micropolar elastic materials undergoing infinitesimal deformation, followed by the principle of minimum potential energy that yields the governing equations to be solved. The corresponding variational form is recovered by the standard procedure while we adopt the Taylor-Hood finite element space for the displacement and microrotation fields that satisfies the Ladyzhenskaya-Babuška-Brezzi (LBB) stability condition (Chapelle and Bathe, 1993;Babuška and Narasimhan, 1997;Sun et al., 2013Sun et al., , 2017Sun, 2015;Na and Sun, 2017). The length-scale-insensitive phase field formulation is then extended to the plane micropolar elasticity, which enables one to simulate cohesive fracture for the higher order continua. We also describe the solution strategy for the micropolar phase field model based on the operator-split, which may potentially be more efficient compared to the monolithic solution scheme.

Mixed Finite Element for Plane Micropolar Elasticity
We first present the implementation of the micropolar elasticity problem.

Problem Statement
Consider a micropolar elastic solid under plane strain condition that occupies a domain B ∈ R 2 with the boundary ∂B. Unlike the Cauchy continua, each material point of a micropolar continuum P may experience both translational in-plane displacement u = [u 1 u 2 ] T and out-of-plane microrotation θ 3 (Eringen, 1966(Eringen, , 2012Ehlers and Volk, 1998;Atroshchenko and Bordas, 2015). In this case, microrotation represents the local rotation of the material point itself, which is independent of the spin component of the velocity gradient. The micropolar strainε and microcurvatureφ are defined as follows:ε where 3 E = 3 E ijk is the Levi-Civita permutation tensor. In this 2D setting, the first and second indices (i and j) of the Levi-Civita tensor take values 1 or 2, while its third index is fixed as k = 3. The nonsymmetric micropolar strain tensor can further be decomposed into symmetric and skew-symmetric parts (i.e.,ε =ε sym +ε skew ) Volume 17, Issue 6, 2019 The definitions of micropolar strainε =ε sym +ε skew and microcurvatureφ can easily be transcribed into Python code by using mathematical operators in the UFL of FEniCS as follows in Listing 1. Compared to Eqs. (1)-(4), one can easily note that the variables u and theta in Listing 1 correspond to the in-plane displacement u and out-of-plane microrotation θ 3 , while .dx(i) denotes the spatial derivative ∂/∂x i . Also, we use as_vector() or as_tensor(), which are the operators for scalar to vector/tensor conversion, in order to specify the return value types of the functions (Alnaes, 2012;Alnaes et al., 2014). We then consider the following linear micropolar elastic stored energy density ψ e (ε,φ): We consider the cases where both elastic moduli C and D are positive definite such that it is necessary to put work into the elastic body to deform it from an unloaded equilibrium configuration. In addition, we consider that the material is isotropic such that elastic moduli C and D can be defined as follows: where λ, µ, κ, and γ are the material constants, that are related to the following engineering material parameters (Diegele et al., 2004;Li and Lee, 2009;Atroshchenko and Bordas, 2015): For a two-dimensional linear elastic isotropic body composed of micropolar continua, at least four elasticity parameters are needed to characterize the elastic constitutive responses. In our implementation, we use the material parameters G, ν, l, and N as inputs. The characteristic length l quantifies the range of couple stress, which implies the size-dependent elastic behavior related to the high-order kinematics. On the other hand, the coupling number N indicates the degree of micropolarity of the material. For example, if N = 0, the micropolar elasticity reduces to the classical elasticity, while N = 1 corresponds to the couple-stress theory (Mindlin, 1962;Mindlin and Tiersten, 1962;McGregor and Wheel, 2014). On the basis of the strain decomposition in Eqs. (3) and (4), our previous study (Suh et al., 2020) decomposed the strain energy density in Eq. (5) into three different parts: (1) the Boltzmann part ψ B e (ε sym ), (2) the microcontinuum coupling part ψ C e ε skew , and (3) the pure microrotational part ψ R e φ . Since we will further incorporate the phase field model to simulate cohesive fracture in micropolar material in Section 2.2 this energy-split strategy enables us to provide a platform that is capable of exploring the effects of distinct degradation on different energy-conjugated pairs. The strain-energy-splitting reads where the partitioned energy densities can be written as follows: ψ C e ε skew = 1 2 κ ε skew :ε skew , The strain energy densities in Eqs. (9)-(11) are used to define force stressσ and couple stressm R , by taking the derivatives with respect to the micropolar strain and microcurvature, respectivelȳ In Eq. (12), we can also decompose force stress into the Boltzmann partσ B and the microcontinuum coupling part σ C , such that σ B ,ε sym , σ C ,ε skew , and m R ,φ become distinct energy-conjugated pairs. The partitioned strain energy densities can thus be recovered as follows: Volume 17, Issue 6, 2019 Similar to Listing 1 our FEniCS code transcribes the force and couple stresses defined in Eqs. (12) and (13)  In Listing 2 tr denotes the trace of a given tensor; Identity(2) corresponds to the second-order identity tensor I for this two-dimensional setting; lamda, mu, kappa, and gamma are the material parameters, which correspond to λ, µ, κ, and γ, respectively. On the basis of the fundamental lemma of calculus of variations, seeking the stationary point of the Lagrange functional yields the Euler-Lagrange equation to be solved. Therefore, assuming no body forces or body couples, and by taking the first variations of Eq. (5) with respect to the field variables (u and θ 3 ), we finally obtain the following set of coupled field equations (Neff, 2006;Li and Lee, 2009;Eringen, 2012;Atroshchenko and Bordas, 2015): We are now equipped with all the necessary ingredients for finite element analysis for the plane micropolar elasticity. Starting from Eqs. (15) and (16), the mixed finite element formulation and the corresponding FEniCS code blocks will be presented in Section 2.1.2.

Mixed Finite Element Formulation
The variational formulation starts with specifying the boundary conditions appropriately. Since the domain boundary ∂B = ∂B u ∪ ∂B tσ = ∂B θ ∪ ∂B tm is composed of Dirichlet (displacement ∂B u and microrotation ∂B θ ) and Neumann boundaries (traction ∂B tσ and moment ∂B tm ), we can specify the boundary conditions as follows: where n is the outward unit normal on the boundary surface andû,θ 3 ,t σ , andt m are the prescribed displacement, microrotation, traction, and moment, respectively. For the purposes of the finite element implementation, we define the trial spaces for the solution variables V u and V θ as, where H α denotes the Sobolev space of order α. The trial spaces in Eqs. (18) and (19) indicate that we employ the Taylor-Hood finite element space to ensure numerical stability (Bathe, 2001;Brezzi and Fortin, 2012). Similarly, the corresponding test function spaces are defined as follows: Finally, the standard weighted residual procedure yields the weak statement for Eqs. (15) and (16), which is to find where Eq. (22) is the weak statement of the balance of linear momentum and Eq. (23) is the weak statement of the balance of angular momentum, respectively. Equipped with Listings 1 and 2 our FEniCS implementation follows the similar procedure described in Eqs. (17)-(23). The corresponding code block starts with defining the finite element function space as illustrated in Listing 3.
Listing 3: Definition of the finite element function space by using mixed element Here, we use VectorElement and FiniteElement to specify the elements. The first argument CG denotes the element type, implying that we adopt the standard Lagrange elements. The second argument mesh.ufl_cell() indicates the unit cell of the discretized domain mesh and the third argument specifies the degree of the element. Then, the operator * constructs the mixed finite element (i.e., the Taylor-Hood element in this case), while the corresponding function space can be defined by using FunctionSpace. On the basis of the function space V that we defined, we specify the trial and test functions [cf. Eqs. (18)-(21)] as follows in Listing 4.

Listing 4: Definition of trial and test functions
Finally, assuming that the problem domain boundary is moment-free, we define the variational form [cf. Eqs. (22) and (23)] by using the predefined trial and test functions in Listing 4. As FEniCS is capable of defining the coupled partial differential equations into one compound system, we directly define the compound variational form as follows in Listing 5.
Listing 5: Plane micropolar elasticity in variational form Here, a is a direct transcription of the bilinear form where the expression inner followed by * dx denotes the inner product, and L is the linear form in which * ds(1) indicates the surface integration over the predefined Neumann boundary. As one may have noted, the close correspondence between the mathematical formulas in Eqs.
(1)-(23) and Listings 1-5 highlights the key strength of FEniCS; the capability of directly translating the mathematical model into the similar Python code (Logg et al., 2012;Alnaes et al., 2015;Langtangen et al., 2016). Once the variational form and Dirichlet boundary conditions (e.g., BC in Listing 6) are properly defined, the only remaining step is to ask FEniCS to compute the solution. This can also be done within few lines of code as follows in Listing 6.

Extension to Micropolar Phase Field Fracture
Using the FEniCS implementation of a two-field mixed FEM for plane micropolar elasticity as our starting point, we now describe the procedure to extend this code to simulate the cohesive fracture in micropolar continua in a phase field framework. We design a code that requires minimal intervention to extend the model, while introducing a modular design that makes it easy to modify the code with different driving force and degradation function for the phase field fracture models.

Phase Field Approximation for Damaged Micropolar Continua
Phase field fracture models employ an implicit density function Γ d (d, ∇d) to indicate the locations of evolving sharp cracks Γ represented by regularized domain (Bourdin et al., 2008(Bourdin et al., , 2012Miehe et al., 2010a) where d ∈ [0, 1] is the phase field or damage variable that varies from d = 0 in undamaged regions to d = 1 in completely damaged regions. In this study, we consider the following general form of the surface density functional (Lorentz et al., 2012;Geelen et al., 2019): where l c is the regularization length that controls the size of the diffusive zone and w(d) is the local dissipation function, which governs the shape of the regularized profile (Clayton and Knap, 2011;Mesgarnejad et al., 2015;Wu and Nguyen, 2018;Wu et al., 2020). Although the quadratic local dissipation w(d) = d 2 has become the standard in modeling brittle fracture, previous studies have pointed out that the linear local dissipation w(d) = d allows to recover the pure elastic phase until the stored energy reaches a certain threshold, which is independent to the regularization length l c (Lorentz et al., 2012;Wu and Nguyen, 2018;Bleyer and Alessi, 2018;Geelen et al., 2019). In order to take advantage of the regularization length insensitive cohesive response, our previous work (Suh et al., 2020) adopted the linear local dissipation function so that Eq. (25) can be simplified into Phase field approximation in Eq. (24) indicates that evolution of the damage variable within a body B corresponds to the generation of new crack surfaces. Therefore, the total potential energy of a damaged micropolar material Ψ is composed of the elastic strain energy [Eq. (5)], and the surface energy that contributes to the crack growth where G c is the critical energy release rate and ψ bulk is the degrading bulk energy of the material. Here, in order to account for a tension-compression asymmetry, we decompose the Boltzmann energy into positive and negative parts (i.e., ψ B e = ψ B+ where • ± = (• ± |•|)/2 is the Macaulay bracket operator,ε sym a , and n a denote the principal Boltzmann strain and its corresponding direction, respectively. In Eq. (27), g i (d) denote the distinct degradation for the corresponding fictitious undamaged energy densities ψ i e (i = B, C, R), i.e., where we assume that the strain energy density parts can either be degraded (i ∈ D) or remain completely undamaged (i ∈ U). In order to replicate the cohesive fracture response, we adopt a special degradation function g(d) that has an associated upper bound on the regularization length scale (Lorentz et al., 2012;Geelen et al., 2019), i.e., where m ≥ 1 is a constant, p ≥ 1 is a shape parameter that affects the material response by controlling the size of the fracture process zone, and ψ crit is the threshold energy density that restricts the crack growth to initiate above this value. From our previous study (Suh et al., 2020) we require that: It is noted that our choice of degradation function g(d) with a specific m parameter plays an important role in phase field approximation of cohesive fracture. However, we omit the justification for these choices for brevity. Interested readers are referred to Lorentz et al. (2012), Geelen et al. (2019), and Suh et al. (2020) for further information.
Up to this point, we prepared the essential parts to extend our FEniCS code presented in Section 2.1 into a phase field model for cohesive fracture in micropolar material. Here, in addition to the micropolar material parameters G, ν, l, and N , we require four more input parameters: G c , l c , ψ crit , and p. As we already defined the stress and strain measures in Listings 1 and 2, we additionally implement the essential components for the extension as follows in Listing 7. Observe that the energy density parts defined in Listing 7 directly inherit the predefined functions described in Listing 1. In this code block, we aim to specify the bulk energy subjected to the degradation (i.e., i∈D ψ i e ) through the if statement with a predefined variable degradation, which is equivalent to set D in Eq. (29). We also define the degradation function g(d) and its derivative [i.e., g_d(d) and g_d_prime(d)], where m and p correspond to the parameters m and p in Eq. (30).
By assuming no body forces or body couples, we finally arrive at the set of Euler-Lagrange equations to be solved, which can be obtained by taking the first variations of the energy functional Ψ with respect to three field variables (u, θ 3 , and d) Following the treatment used in (Miehe et al., 2015;Bryant and Sun, 2018), we ensure the crack irreversibility by introducing a history function H. In particular, we define H as a pseudo-temporal maximum of the degrading strain energy density and at the same time as a variable that restricts the crack growth to initiate above a threshold ψ crit to approximate the cohesive response. On the basis of this treatment, Eq. (34) becomes where the definition of the history variable H can be transcribed into a simple Python code block as follows in Listing 8. As illustrated in Listing 8, function H is designed to return a pseudo-temporal maximum bulk energy by comparing the current bulk energy with the previous history variable through the UFL operator conditional (Alnaes, 2012;Alnaes et al., 2014). We now collected all the additional components to extend our FEniCS code for plane micropolar elasticity presented in Section 2.1. Also starting from the Euler-Lagrange equations in Eqs. (32), (33), and (35), the variational formulation and the corresponding FEniCS implementation will be presented in Section 2.2.2.

Variational Formulation
Recall Section 2.1.2 that the domain boundary can be decomposed into Dirichlet (∂B u and ∂B θ ) and Neumann (∂B tσ and ∂B tm ) boundaries, we thus prescribe the boundary conditions as follows: The initial conditions are imposed as u = u 0 , θ 3 = θ 3 | 0 .
Since we now have three field variables (i.e., u, θ 3 , and d), total three trial spaces should be defined, i.e., The corresponding test function spaces read Equations (38)-(43) imply that we employ the Taylor-Hood finite element space for u and θ 3 while the phase field d is discretized via linear finite elements. By adopting similar process employed to obtain Eqs. (22) and (23), we get the weak statement for the Euler-Lagrange equations in Eqs. (32), (33), and (35); that is, find {u, where G u → R is the weak statement of the balance of linear momentum G θ → R is the weak statement of the balance of angular momentum and G d → R is the weak statement of the damage evolution equation: The corresponding FEniCS transcription again starts with defining the finite element function spaces. Since we have already defined the Taylor-Hood finite element function space for the displacement and microrotation fields in Listing 3, our new task is to define the function spaces for the phase field and the history variable as follows in Listing 9.
Listing 9: Definition of the function spaces for the phase field and the history variable 1 # Define function spaces 2 W = FunctionSpace(mesh, 'CG', 1) # phase field 3 WW = FunctionSpace(mesh, 'DG', 0) # history variable Here, CG and DG in the second argument of FunctionSpace denote the standard Lagrange element and discontinuous Lagrange element families, respectively. Once we define the function space appropriately, the next step is to define trial and test functions. Since the damage evolution equation in Eq. (47)  As we later employ a solution scheme based on the operator splitting (which will be described in Section 2.2.3), we intentionally define two variational problems that are decoupled. Unlike the linear variational problem in Listing 6, each system requires us to define the Jacobian matrices to specify the problem through NonlinearVariational Problem. However, one of the strengths of FEniCS is that it is capable of symbolically deriving the Jacobian matrix by using the operator derivative as illustrated in Listing 11, so there is no need to manually derive the Jacobian matrix explicitly.

Operator-Split Solution Scheme
We now describe an efficient solution procedure and its FEniCS implementation based on the operator-split solution scheme, which is also known as the staggered scheme. As previously hinted by Listing 11, the main idea in this approach is to algorithmically decouple the system and successively update the field variables (Miehe et al., 2010a;Na and Sun, 2018;Bryant and Sun, 2018;Geelen et al., 2019). In the first part of this operator-split setting, we update the phase field while the displacement and microrotation fields are held fixed. Since the damage evolution equation in Eq. (35) is nonlinear, the phase field is obtained iteratively once the nonlinear solver converges within a predefined tolerance. The second part then advances the displacement and microrotation fields with a new damage field being fixed. Compared to the monolithic approach, this operator splitting may potentially be more efficient and robust as pointed out in Miehe et al. (2010a).
As we described above, the staggered scheme requires us to adopt two different solvers for each system. Thus, a necessary preprocessing step for the solution scheme in FEniCS is to specify the solvers as follows in Listing 12. We employ the Newton-Raphson method to advance the displacement and microrotation fields, while the Scalable Nonlinear Equations Solvers (SNES) from the open-source toolkit PETSc is adopted to solve the damage evolution equation since the phase field possesses lower and upper bounds (Balay et al., 2001(Balay et al., , 2019. One can also tune the solvers by specifying the solver parameters as illustrated in Listing 12. Now we have arrived at our final destination: the FEniCS implementation of the operator-split solution scheme. Unlike Listing 6, this solution scheme should be programmed explicitly. As illustrated in Listing 13, the most outer loop corresponds to the increasing sequence of discrete times from t_i to t_f. At each time step, we prescribe the displacement-driven incremental load on the specified region, followed by the staggered loop. The staggered iteration solves the phase field equation through solver_d.solve() to update the phase field, and then solves the balance equations to advance the displacement and microrotation by solver_ut.solve(). This process is repeated until the residual reaches the predefined tolerance staggered_tol or the number of iteration reaches the specified staggered_maxiter.
Listing 13: Operator-split solution scheme 1 # Staggered scheme 2 t = t_i 3 while t <= t_f: In this section, we derived a phase field model for cohesive fracture in micropolar continua step by step and provided their counterparting FEniCS code blocks at each step. Since we only described the partitioned code blocks throughout Section 2.2 we complete this section by providing a flow chart that summarizes the entire code structure. As illustrated in Fig. 1, the very first step to solve a boundary value problem is to specify all the input parameters (i.e., material properties G, ν, l, N , G c , l c , ψ crit , and p, as well as solver parameters) and importing the mesh. Having defined the boundary conditions on domain boundaries properly, the program defines the internal variables as we previously described in Listings 1, 2, 7, and 8. Once the variational form is defined, our implementation then solves the system of equations by adopting operator-split approach discussed in Section 2.2.3.

NUMERICAL EXAMPLES
The purpose of this section is to verify the implementation with available benchmarks and demonstrate the capability of the FEniCS implementation of a phase field fracture model for micropolar continua. Our first numerical example examines the classical problem of stress concentration around the hole. As a verification exercise, we compute the stress concentration factors and the stress field inside the domain with different coupling numbers N and compare to the analytical solutions. The second example aims to showcase the performance of the code by studying crack propagation in double notched micropolar elastic specimen subjected to combined tensile and shear loads. All the

Verification Exercise: A Circular Hole in an Infinite Plate
If an infinite plate (composed of nonpolar continua) with a circular hole is subjected to remote uniaxial tension, the tangential stress reaches a value three times the applied stress at two points on the edge of the hole perpendicular to the loading axis (Green, 1940). This factor is often referred to as the stress concentration factor. Ariman (1967) has extended this stress analysis to micropolar elasticity material and find that this stress concentration factor may vary, depending on the micropolar effect for micropolar continua. Our first numerical example serves as a benchmark problem for verifying the numerical implementation of the mixed finite element code for the micropolar elastic material introduced in Section 2. In particular, we would like to verify whether we can recover the stress concentration factors obtained from analytical solutions in the literature.
We consider an rectangular plate under pure tension σ 0 that has a circular hole with diameter D at the center. By symmetry, the problem domain is modeled as a finite quarter-plate of size 5D × 5D instead, which is assumed to be sufficiently large enough to neglect any possible boundary effects (Fig. 2). As for the boundary conditions, left and bottom boundaries are supported by rollers with fixed microrotation, while we prescribe a uniform tension σ 0 on the right boundary. All other boundaries are maintained traction-and moment-free during the simulation. Also, we choose the following material parameters for this example, shear modulus G = 50 GPa, Poisson's ratio ν = 0.3, and bending characteristic length l = D. Within this problem domain, we investigate the effect of micropolar effect on the stress concentration factors (SCF) (i.e., the ratio between the highest force stress to the applied stress σ 0 ) and the distribution of the stress fields by considering different coupling numbers N = 0.1, 0.25, 0.5, 0.75, and 0.9. Figure 3 shows the exemplary distribution of normalized stress fields within the domain when the coupling number is set to be N = 0.9. As expected, the maximum stress occurs on the edge of the hole perpendicular to the  (Anderson, 1991;Yang et al., 2008;Eringen, 2012;Atroshchenko and Bordas, 2015). Table 2 summarizes the numerically computed stress concentration factor (SCF FEM ) compared to the analytical solution (SCF analytical ) given in Eringen (2012). Here, we define the relative error (e) as follows: The overall results in Table 2 demonstrate that the stress concentration factor for the plate with a circular hole can successfully be reproduced by our FEniCS implementation within 2.7% of relative error. The numerical and analytical results also indicate that the maximum force stress becomes smaller as the coupling number N increases, while the (1) and force stress in Eq. (12) imply that the normal stresses are equivalent to those in the classical approach, whereas the shear stresses depends on the microrotation. We therefore compare the analytically and numerically obtained shear stress distributions along OA in Fig. 2, where the shear stress components in polar coordinates (i.e.,σ rθ andσ θr ) are extreme. Figure 4 shows the distribution of the normalized shear stress components along OA (i.e., θ = 45 deg in polar coordinates), where the solid curves with symbols indicate the numerical results obtained from our FEniCS implementation, while the transparent curves denote the analytical solutions from Eringen (2012). Since we set the bending characteristic length l = D to be small compared to the size of the domain (5D × 5D), all the results tend to asymptotically converge toward the value of -0.5 as r increases regardless of the degree of micropolarity. Furthermore, in the region where r < l the degree of skewness in force stress tensor (i.e., the difference between shear stress components) tend to increase as N increases, which implies that the coupling number N directly controls the degree of asymmetry of the model. The numerical results in good agreement with the analytical solution verify the correctness of the FEniCS implementation of mixed finite element method for micropolar elasticity that we developed in this study.

Double-Edge-Notched Test
Since the detailed investigation on the micropolar phase field model (including regularization length insensitive material response, partial degradation, and micropolarity effects on the crack pattern) can be found in our previous study (Suh et al., 2020), this numerical example showcases the capability of our FEniCS implementation for capturing crack nucleation, propagation, and coalescence in micropolar material. Consider a 100 mm wide and 100 mm long square plate with two 25 mm long edge notches as illustrated in Fig. 5. The problem domain is subjected to combined tensile and shear loads, by prescribing the displacementū along the top boundary at an angle of 45 deg with ∆u 1 = ∆u 2 = 5.0 × 10 −4 mm, while bottom boundary is held fixed. By assuming that all the energy density parts can be degraded (i.e., D = {B, C, R} and U = ∅), the material parameters for this example are chosen as follows: shear modulus G = 12.5 GPa, Poisson's ratio ν = 0.2, bending characteristic length l = 30.0 mm, coupling number N = 0.5, critical energy release rate G c = 0.1 N/mm, threshold energy density ψ crit = 0.1 kJ/m, and the shape parameter p = 10. Figure 6 illustrates the evolution of the phase field at several load increments, while Fig. 7 shows the global response of the structure during the simulation. Recall Eqs. (26) and (30)  In this case, unlike standard phase field model for brittle fracture that adopts w(d) = d 2 and g(d) = (1−d) 2 , the crack is expected to propagate even without complete separation, since our model is capable of considerably elongating the length of the fracture process zone by controlling the parameter p. As expected, we observe the partially damaged region (i.e., 0 < d < 1) in Fig. 6, where complete failure may eventually take place herein. The crack starts to nucleate at notch tips as the stored energy density exceeds the predefined threshold ψ crit . Then, the phase field starts to propagate atū 1 =ū 2 ≈ 0.01 mm, where the material enters the softening deformation stage (Fig. 7). Compared to the results in Bocca et al. (1991) and Wang and Sun (2017), where neighboring two cracks do not coalesce, in this example, two cracks kink toward the adjacent notch and eventually coalescence with each other. As Yavari et al. (2001) pointed out, crack propagation in micropolar materials requires additional effort on breaking the microrotational bond. Since cracks generally tend to propagate in the direction that maximizes the dissipation (i.e., minimizing the additional bond-breaking effort in this case), the observed crack trajectories are consistent with this interpretation. The results underscore that introducing size-dependent constitutive model itself may affect the crack trajectory, while we also highlight that our FEniCS codes are capable of capturing the micropolar effects on both elastic and damaged regimes.

CONCLUSION
We provide a detailed account of an implementation of phase field fracture models for micropolar materials via FEniCS. By revisiting the theory established in Suh et al. (2020) this study presented a brief step-by-step overview of the mathematical model and at the same time provided a detailed description on the design and implementation of the computer code that runs the mathematical model. We first present the procedure to implement the mixed FEM for plane micropolar elasticity. We then extend this platform into a phase field code that simulates cohesive fracture therein by explicitly programming the staggered solution scheme. We also provide numerical examples to verify the correctness and showcase the capability of our program. The key strength of this implementation is the minimal approach enabled by the FEniCS libraries and the corresponding Application Programming Interfaces (API), which enable us to write Python code with expressions easily understood by end users without comprising significantly on efficiency. By making the source code available in a repository, this work will enable third-part verification and validation and therefore promote transparency and collaborations. The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the sponsors, including the Army Research Laboratory or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Government purposes notwithstanding any copyright notation herein.