Multiscale dynamic fracture analysis of composite materials using adaptive microstructure modeling

In this study, a computational framework is proposed to investigate multiscale dynamic fracture phenomena in materials with microstructures. The micro‐ and macro‐scales of a composite material are integrated by introducing an adaptive microstructure representation. Then, the far and local fields are simultaneously computed using the equation of motion, which satisfies the boundary conditions between the two fields. Cohesive surface elements are dynamically inserted where and when needed, and the Park‐Paulino‐Roesler cohesive model is employed to approximate nonlinear fracture processes in a local field. A topology‐based data structure is utilized to efficiently handle adjacency information during mesh modification events. The efficiency and validity of the proposed computational framework are demonstrated by checking the energy balances and comparing the results of the proposed computation with direct computations. Furthermore, the effects of microstructural properties, such as interfacial bonding strength and unit cell arrangement, on the dynamic fracture behavior are investigated. The computational results demonstrate that local crack patterns depend on the combination of microstructural properties such as unit cell arrangement and interfacial bonding strength; therefore, the microstructure of a material should be carefully considered for dynamic cohesive fracture investigations.

objective of the present study is to develop an efficient and effective computational framework, which can account for the nonlinear dynamic fracture behavior associated with spontaneous multiple microcrack initiation and branching in conjunction with the microstructure of a material.
To investigate a material with microstructure, various theoretical and computational methods have been developed, including multiscale approaches. 4 An average macroscopic constitutive relation was obtained from microstructural constituents using the concept of mathematical and computational homogenization. [5][6][7][8][9][10] However, loss of material stability due to cracks or localization results in strain softening, and thus the standard homogenization rules, for example, periodic boundary conditions and volume average of stress, are not applicable for general fracture problems. 7,11 To account for material failure with respect to microstructure, for example, the loss of material stability, various multiscale computational techniques have been proposed. 12 One main category of multiscale methods is named as semiconcurrent method. In this approach, two scales (eg, coarse and fine) are generally introduced, and fine and coarse scale problems were solved using various computational methods, for example, finite element method, generalized/extended finite element method, discrete element method, molecular dynamics, and virtual element method. [13][14][15][16][17] When both fine and coarse scale problems are solved with finite element method, this approach is also called as FE 2 method. 18,19 To provide the consistent fracture energy between two scales, for example, regularization parameters were transferred from fine scale to coarse scale. [20][21][22] Furthermore, to satisfy the compatibility and equilibrium conditions between two scales, several approaches were proposed in conjunction with domain decomposition and bridging zone, [23][24][25][26] named as concurrent method.
In the multiscale computational methods, one should note that a smooth transition between micro-and macro-scales is essential, especially for wave propagations. This is because when a wave travels between microand macro-models, spurious wave reflections occur, which can affect the crack speed, propagation direction, and patterns. 27,28 Furthermore, dynamic fracture is generally associated with spontaneous crack initiation and branching, which leads to complex crack patterns, and thus microscale information (eg, length scale, boundary conditions, effective displacement jump, traction-separation relations, and others) should be carefully transferred to the macroscale.
In continuum mechanics, nonlinear material fracture or failure at micro-and/or macro-scales is generally represented by either continuum or discrete approaches. For the continuum approaches, strain localization and material softening are introduced within narrow bands in conjunction with intrinsic length-scale measures in the constitutive theory. [29][30][31] A recent phase-field model of fracture can be classified as a special case of a gradient plasticity model. [32][33][34] In the case of discrete approaches, the concept of the cohesive zone model is utilized to approximate the nonlinear fracture process zone, while discontinuities due to fracture are explicitly described within a continuum domain using cohesive surface elements, [35][36][37][38][39] which is the model of choice in the present study. Alternately, discontinuities on a continuum domain can be introduced using enrichment functions in conjunction with the partition of unity, named as GFEM/XFEM. [40][41][42] Further in-depth reviews on the damage mechanics and the cohesive zone model can be found in the literature. 32,35,43,44 In this study, to consistently handle the boundary between micro-and macro-models for analyzing the dynamic cohesive fracture of a composite material, a multiscale computational framework is proposed using an adaptive microstructure representation. Based on the macroscale information, a local field is solved, and the solution of the local field is directly embedded in a macromodel. Then, the equation of motion is simultaneously solved at the micro-and macro-scales, owing to which the proposed framework is not affected by the boundary conditions between micro-and macro-models. Fracture in the micromodel is explicitly described by adaptively inserting cohesive surface elements in a local domain, and the traction-separation relation is defined using a potential-based model. The proposed computational framework is verified by checking the energy balance and comparing the results of adaptive modeling and direct computations. Furthermore, the computational results demonstrate that the microstructure of a material can significantly affect its dynamic cohesive fracture characteristics, such as the dissipated energy, crack velocity, and crack patterns.
The rest of the article is organized as follows. In Section 2, a multiscale computational framework is presented by integrating micro-and macro-scales in conjunction with an adaptive microstructure modeling scheme. Computational implementation is described in Section 3. In Section 4, the proposed computational framework is verified by checking the energy balance and comparison between adaptive modeling and direct computations. Section 5 demonstrates the effect of microstructure on dynamic cohesive fracture characteristics, such as the dissipated energy, crack velocity, and crack patterns. Finally, the key findings of the present study are summarized in Section 6.

INTEGRATION OF MICRO-AND MACRO-SCALES
To integrate the macro-and micro-scales of a composite material, an adaptive microstructure modeling scheme is proposed. The microstructure of a material is explicitly considered to accurately capture localization at potential hot spots (eg, initial defects or cracks), and effective material properties are employed to represent macroscopic material behavior. When a crack propagates, the microstructure of a material is adaptively embedded and local material properties are consistently updated to bridge the micro-and macro-scales.

Transition from a far field to a local field
A homogenized material at a far field is adaptively transformed to a material with a microstructure according to the crack propagation. A coarse discretization is used for the far field and is adaptively refined to represent the local field. In adaptive local field representation, new nodes and elements are created during dynamic fracture analysis. Thus, material properties and physical field quantities are consistently updated for new nodes and elements. Material properties are simply assigned to each new element for a given microstructure because the boundaries of the particles correspond to the boundaries of finite elements in this study. Physical quantities, that is, velocities and displacements of new nodes, are updated when a homogenized material (ie, far field) is transformed to a material with microstructure (ie, local scale). First, the velocities of new nodes are evaluated by interpolating from existing nodal velocities. Lagrange-basis shape functions are employed for the interpolants and thus the velocity vector at a given point (x n , y n ) can be written as v(x n , y n ) = Nv, (1) where N(r n , s n ) is a shape functional matrix and v is a vector of nodal velocities. The position of a new node in natural coordinates (r n , s n ) is obtained from the inverse mapping of a coordinate transformation (T −1 t ), that is, Transformation from the natural coordinates to the physical coordinates is simply defined using the Lagrange-basis shape functions of an element. Then, the inverse mapping of T t is performed using Newton's algorithm.
Displacements of new nodes are approximated by introducing local analysis of a unit cell, as shown in Figure 1. A local model of the unit cell is first constructed by adaptively refining a coarse mesh, as described in Section 2.3. Next, displacement along the edges of a unit cell is evaluated by interpolating the nodal displacements with Lagrange-basis shape functions. Then, local finite element analysis is performed with the displacement boundary conditions along the edges of the unit cell. The solution of the local model is utilized to update the nodal displacement within a unit cell.
Alternatively, the nodal displacements of new nodes can be interpolated from existing nodes using Lagrange-basis shape functions, which is the same approach as nodal velocity evaluation. One should note that the use of Lagrange-basis shape functions leads to strain and kinetic energy conservations when adaptive mesh refinement is performed for homogeneous materials. 45 When a homogenized material transforms into a composite material during mesh refinement, the interpolated displacement can result in relatively large errors in the local deformation and stress fields, which will be discussed in Section 4.1.
In summary, the physical quantities, such as velocities and displacements, in a local field are updated in adaptive microstructural modeling. Subsequently, the local and far fields are solved at the same time using the equation of motion, which leads to a consistent evaluation of the displacement, velocity, and acceleration in a time integration scheme, as discussed in Section 3.1. It is noted that the approximation of acceleration for new nodes is not needed because acceleration is evaluated using the equation of motion after adaptive microstructural representation.

Far field approximation
To approximate the elastic behavior of composite materials in a far field, a simple micro-mechanical model, that is, the Mori-Tanaka method, 46 is employed. The effective material properties of a far field are evaluated using the properties of the matrix and particle at a local field with a given volume fraction (f ). The effective bulk modulus (K) and effective shear modulus ( ) are given by and where K m and K p are the bulk moduli of the matrix and particle, respectively, and m and p are the shear moduli of the matrix and particle, respectively. Additionally, the densities of the matrix and particle ( p and m ) are transformed into an effective density for the far field using the rule of mixtures as follows which leads to mass conservation between the local and far fields.

Local field representation
To describe the microstructure of a composite material in a local field, a finite element mesh is adaptively refined while cracks propagate. For mesh refinement, the element size should be smaller than a certain length scale, which is associated with the size and shape of the microstructure. In this study, a circular particle is employed in the material microstructure. Furthermore, to describe circular particles in a local refined mesh, local refined nodes are relocated to fit with particle boundaries. Figure 2 illustrates a procedure corresponding to local field representation. A relatively coarse 4k mesh is employed for far field representation, as shown in Figure 2A. A 4k mesh consists of right isosceles triangles, which include two legs and a hypotenuse. Then, the 4k mesh is refined by splitting the hypotenuses, which are termed as edge split operator. 45 The edge split operation can be recursively conducted because it does not change the topology of the 4k mesh while decreasing the element size ( Figure 2B). Next, within a refined domain, a circular particle is represented by relocating local nodes, as shown in Figure 2C. Finally, a nodal perturbation is introduced to reduce mesh-bias in the 4k structured mesh 47 and a Laplacian operator is utilized to improve the mesh quality. Figure 2D illustrates a local finite element mesh with a circular particle.
To explicitly and accurately describe the particle shape in a refined mesh (eg, Figure 2C), adequate edges, which correspond to a particle boundary, are identified using the following process (see Figure 3). First, a starting point is arbitrarily selected on the particle boundary. Then, a vertex (or corner node) of a triangular element, which is nearest to the starting F I G U R E 2 Adaptive local field representation: A, Initial far field, B, adaptive mesh refinement, C, particle boundary representation, and, D, nodal perturbation and Laplace operator point, is identified and the identified vertex is relocated to the starting point. In Figure 3A, the particle boundary is illustrated as a dashed curve, while the starting and relocating position vectors are denoted as v s andv, respectively. Next, from the starting point, the adjacent edges and their vectors are calculated (ie, e i = v i − v s ), where the ith adjacent edge connects the starting point (v s ) and the ith adjacent point (v i ), as shown in Figure 3B. Then, the angles between the edge vector (e i ) and unit tangent vector to the particle boundary (n t ) are evaluated, and the edge, which provides the smallest angle, is selected as an adequate edge. For instance, e 2 corresponds to the adequate edge vector in Figure 3B and thus v 2 is projected and relocated to the particle boundary (see Figure 3C). Then, the relocated position vector (v) becomes a new starting point and an adequate edge vector from a new starting point is identified, as shown in Figure 3D. The adequate edge search procedure is recursively implemented until a closed domain is formed. Local mesh modifications, which include vertex projections to the particle boundary and nodal perturbation, can result in low-quality elements for finite element computation. To improve the mesh quality, a Laplacian operator is utilized. The Laplacian operator of a given vertex (v) is defined as where i is the weight of adjacent vertices (v i ) and n e is the number of the adjacent edges of the vertex (v). In this study, every edge vector has the same weight and thus i is 1/n e . Then, the vertex (v) is relocated to a new vertex (v ′ ), that is, One should note that the Laplacian operator is not employed for vertexes associated with particle boundaries in order to maintain the particle shape.
To further improve the mesh quality, some of the vertexes of the particle boundaries are updated. Specifically, when the two legs of an element are adequate edges of a smooth particle boundary ( Figure 4A), all three nodes of a triangular element are located at the particle boundary and thus the Laplacian operator cannot be employed in this case. Furthermore, the angle between the two legs is large for a smooth boundary, which results in a lower mesh quality ( Figure 4B). To tackle such case, the particle boundary is represented by the hypotenuse of an element and the node, which shares the two legs, is relocated on the basis of the Laplacian operator, as shown in Figure 4C.

COMPUTATIONAL IMPLEMENTATION
The finite element method is formulated based on a reference configuration, that is, the total Lagrangian formulation, while the central difference method is employed to solve transient problems. A potential-based cohesive zone model 48 is utilized to approximate the nonlinear fracture process zone. Discontinuities like cracks in a continuum domain are adaptively introduced by inserting surface elements in conjunction with a topology-based data structure 39,49,50 which leads to extrinsic cohesive zone modeling. 47,51

Finite element formulation
The weak form of mechanical problems is obtained from the principle of virtual work, which is given as where E and S denote the Green-Lagrange strain tensor and the second Piola-Kirchhoff stress tensor, respectively, and δu represents virtual displacement. The external traction (T ext ) acts on the boundary surface (Γ 0 ), and the traction vector (T fra ) on the fracture surface (Γ fra ) consists of cohesive traction (T coh ) and contact force contribution (T con ).

Time integration procedure for adaptive microstructure modeling
A central difference method, which is a special case of the Newmark-method, is utilized. The time integration procedure for adaptive microstructure representation is summarized in Algorithm 1. After applying boundary conditions, the velocity (u n+ 1 /2 ) at time step (n+ 1 / 2 ) and the displacement (u n+1 ) at time step (n + 1) are evaluated for the time increment of Δt where the superimposed dot denotes the time derivative. Based on the crack initiation criterion in conjunction with the stress state at time step (n + 1), new fracture surface is adaptively created using surface elements. Around the new fracture surface, a local microstructure is explicitly described with adaptive local field representation. The material properties are updated, and the velocities and displacements of the new nodes are consistently evaluated, as discussed in Section 2. Next, the external force (f ext,n+1 ) and internal force (f int,n+1 ) at time step (n + 1) are calculated. On the fracture surface, the cohesive force (f coh,n+1 ) is obtained from a cohesive traction-separation relationship for the time step (n + 1). Additionally, to prevent from material self-interpenetration, the contact force (f con,n+1 ) is estimated using a simple penalty method 57 in conjunction with a normal gap at time step (n + 1) and the normal gap rate at time step (n + 1 / 2 ). The detailed estimation procedures of f coh,n+1 and f con,n+1 are explained in the following subsection. Finally, the accelerations at time step (n + 1) are evaluated using the equation of motion. The mass matrix (M) is diagonalized by scaling the diagonal terms with respect to the total mass; this is known as the diagonal scaling method. 52 One notes that the equations of motion for the local and far fields are solved at the same time in adaptive microstructure representation.

Traction-separation relationship on fracture surface
Nonlinear material fracture process is represented by the concept of the cohesive zone model. 53,54 The cohesive zone model defines a traction-separation relationship across the fracture surface, which describes progressive material failure. A Park-Paulino-Roesler (PPR) potential-based cohesive model 48 is utilized; it yields a consistent traction-separation relationship for an arbitrary separation path. 35,55 The PPR cohesive model has three physical fracture parameters in each fracture mode, that is, fracture energies ( n , t ), cohesive strengths ( max , max ), and shape parameters ( , ), for the extrinsic cohesive zone model. The potential function of the extrinsic cohesive zone model is defined as The normal and tangential cohesive tractions are obtained from the gradients of the potential and can be respectively expressed as and The characteristic parameters ( n , t , n , and t ) are explicitly determined using fracture parameters, such as the fracture energy, cohesive strength, and shape. 48,56 To avoid material self-interpenetration during dynamic fracture, a simple penalty method is employed with nondimensional penalty parameters. 57 In the penalty method, the normal contact force (f con N ) is evaluated using the normal gap (g N ) and normal gap rate (̇g N ) on the fracture surface as follows where M + and M − are the masses on the upper and lower surfaces, respectively. One notes that p 1 and p 2 are nondimensional contact parameters, which are less sensitive to changes in the material properties and time increments. Based on a previous parametric study, 57 p 1 and p 2 are selected as 0.8 and 0.5, respectively.

Topology-based data structure
A topology-based data structure, called TopS, is employed to retrieve adjacent information in a finite element mesh. 49,50 In the data structure, nodes and elements are explicitly stored in a memory space, while vertices, edges, and facets are described by references to explicit entities. Subsequently, adjacency relationships, for example, node-to-element, edge-to-element, and element-to-element, are efficiently obtained using TopS. Additionally, to maintain a consistent topology-based data structure, a client-server approach is utilized. 47 Based on TopS, a 4k structured mesh is recursively refined using edge split operators for adaptive local microstructure representation. When a surface element is inserted, the nodes are duplicated and local element connectives are updated. A callback function is introduced to obtain access to the duplicated nodes and update nodal quantities.

VERIFICATION
The proposed computational framework is verified by investigating the effects of microstructural representation on the energy balance, elastic response, and fracture behavior. First, kinetic and strain energy differences between local and far fields are quantified for arbitrary displacement and velocity fields. Next, an elastic wave propagation example is studied by comparing the three cases of material representations, that is, a homogenized material, material with a complete microstructure, and material with adaptive microstructure representation. Finally, cohesive fracture computation with the adaptive microstructure representation is compared with the result obtained by direct computation.

Energy balance between local and far fields
Due to material characteristic transitions from a homogenized material to a material with a microstructure, the strain and kinetic energies undergo a change. To quantify energy differences between the local and far fields, the strain and kinetic energies in a local field (ie, material with a microstructure) are compared with the energies in a far field (ie, homogenized material) for given displacement and velocity fields. For a local field, a square domain of 1 mm × 1 mm is selected and a circular particle with a diameter of 0.6 mm is placed in it, which leads to a volume fraction of 0.2827. The elastic modulus and Poisson's ratio of the matrix are 100 GPa and 0.15, respectively, and the density is 2000 kg/m 3 . The elastic modulus, Poisson's ratio, and density of the particle are 150 GPa, 0.2, and 3000 kg/m 3 , respectively. The local microstructure is homogenized for the far field approximation using the Mori-Tanaka method (Equations (3) and (4)) and the rule of mixtures (Equation (5)). For strain-energy evaluation, two arbitrary displacement fields are selected, that is, linear and quadratic fields, as listed in Table 1. Then, the far field is discretized into four constant strain triangular elements, that is, T6, and nodal displacements are assigned from the analytical displacement field of Table 1. For local field evaluation, the domain is discretized into 1024 triangular elements and nodal displacements are obtained from the far field information using either the local analysis of a unit cell or the interpolation with the Lagrange-basis shape functions, as discussed in Section 2.3. The calculated strain energies of the far and local fields are summarized in Table 2. Strain energy difference between the far and local fields is about 2% when local displacement is obtained from the interpolation with the Lagrange-basis shape function, while the difference significantly decreases to 0.5% to 0.7% when the nodal displacement is obtained from the solution of the local field analysis.
Additionally, the horizontal displacements of the far and local fields for the linear case are illustrated in Figure 5. The displacement of the far field ( Figure 5A) is the same as the displacement of the local field obtained from interpolation ( Figure 5B), as expected. However, the displacement of the local field obtained from the solution of the unit cell analysis ( Figure 5C) is slightly different from the far field displacement. Based on the obtained displacement fields, horizontal stress within the domain is evaluated, as shown in Figure 6. Figure 6A exhibits a constant stress field because of the homogenized material and linear displacement field. Figure 6B demonstrates a constant stress field within the particle and matrix, which may not correspond to a realistic stress field for a material with microstructure. When the local displacement is obtained from the solution of the unit cell, the corresponding stress field ( Figure 6C) is more accurate than the constant stress fields ( Figure 6A,B) for a material with microstructure.
For kinetic energy evaluation, linear and quadratic velocity fields are employed by dividing the displacement fields by 0.1 μs; in addition, a consistent mass matrix is used. Similarly, four T6 triangular elements are used for far field discretization, while 1024 elements are used for local field discretization. The calculated kinetic energies are summarized in Table 2 and the kinetic energy difference between the far and local fields is about 2%.

Elastic wave propagation
A cantilever example with microstructure is employed to investigate the effects of adaptive microstructure modeling on elastic wave propagation. The geometry and boundary conditions of the cantilever example are shown in Figure 7. The length and height of the cantilever are 2 and 0.3 m, respectively, while its thickness is 0.1 m. The microstructure of a material is represented by square unit cells. The size of each unit cell is 100 mm by 100 mm and a circular particle is included in each of the unit cells. The diameter of each particle is 60 mm and the volume fraction is 0.2827. Sixty square unit cells are uniformly distributed within the domain. The elastic modulus, Poisson's ratio, and density of the matrix are 20 GPa, 0.15, and 2200 kg/m 3 , respectively. The elastic modulus, Poisson's ratio, and density of the particles are 40 GPa, 0.25, and 2500 kg/m 3 , respectively. Additionally, a sinusoidal transient load is applied at the end of the cantilever, that is, P(t/T) = P 0 sin( t/T), during the natural period of the cantilever (ie, T = 0.8017 second), where t represents time and P 0 is the amplitude of the applied force (ie, P 0 = 1 N). The cantilever example with the microstructure is solved using three approaches, that is, direct computation, homogenized model, and adaptive microstructure model. For direct computation, 60 unit cells are explicitly described in a finite element mesh, and the material properties of the matrix and particles are directly assigned. A magnified view of the finite element mesh around the support of the cantilever is shown in Figure 8A. The element size is about 6 mm and the numbers of elements and nodes are 47 352 and 95 073, respectively. For the homogenized model, the material properties are approximated using the Mori-Tanaka method and rule of mixtures, which can significantly reduce computational cost with a coarse mesh. However, to avoid errors associated with element size in the homogenized model, a similar element size is utilized; the numbers of elements and nodes are 15 360 and 31 089, respectively. The corresponding magnified view of the mesh on the support is shown in Figure 8B. For the adaptive microstructure model, the unit cells are sequentially described from the fixed boundary to the free end in 1 second. A magnified view of the finite element mesh at the initial configuration is shown in Figure 8C and the corresponding number of elements and nodes are 5169 and 10 438, respectively. These numbers eventually increase to the node and element numbers for the direct computation. The displacement at the loading point in the three computational approaches is plotted in Figure 9. The numerical results match well with each other although the result of the homogenized model is slightly different from the results of the other two cases. Additionally, the strain energy density at 0.35 second is illustrated in Figure 10. The overall strain energy distribution is similar in all the cases; a relatively high strain energy is observed at the top and bottom of the support regions. However, the homogenized model cannot capture the local features of the microstructure and thus results in a smooth strain energy distribution. When the microstructure is considered in the computation, a local stress (or strain) redistribution is expected. For example, the matrix exhibits a relatively higher strain energy than the particle because the elastic modulus of the matrix is lower than the elastic modulus of the particle. Such local features can significantly impact the crack propagation behavior. In summary, the global elastic responses of the homogenized material are similar to the responses of the material with microstructure, while a local stress redistribution due to the presence of microstructure can be captured using the adaptive microstructure representation.

Dynamic fracture computation
To verify the adaptive microstructure modeling scheme for dynamic fracture computation, a single-edge notched tension example is utilized. A rectangular domain of 6 mm by 2.8 mm is employed with an initial notch length of 1.55 mm, as with v 0 =, A, 25, B, 80, and C, 250 mm/s and 300 MPa, respectively. The shape parameters in the PPR model are selected as 2, which leads to a linear softening model. To represent the composite material, square unit cells with circular particles are defined within the domain. The size of the unit cell is 400 μm and the radius of the particle is 120 μm, which leads to a particle volume fraction of 0.283. Additionally, a small randomness is introduced for the positions of the particles in the unit cells.
The single edge notched tension example is computed using two approaches, that is, direct computation and adaptive microstructure modeling. For direct computation, 105 unit cells are placed from the initial configuration; a magnified view of the corresponding finite element mesh around the initial notch is shown in Figure 12A. The numbers of elements and nodes are respectively 107 520 and 215 869. In the adaptive microstructural model, a unit cell is placed at the initial notch ( Figure 12B). The numbers of elements and nodes are 2155 and 4392, respectively, at the initial configuration; these numbers are significantly smaller than the numbers for the direct computation. When the crack propagates, the microstructure of the material is adaptively described only around crack tip regions. The computational results demonstrate that the overall crack patterns of adaptive microstructure modeling ( Figure 12D) agree well with the patterns of direct computation ( Figure 12C). The vertical stress fields of direct computation and adaptive microstructure modeling are illustrated in Figures 13 and 14, respectively. A high stress concentration around the crack tip region is observed in both computational approaches, while some features of the microstructure in regions far from the crack tips (see Figure 13) are averaged in adaptive microstructural modeling (see Figure 14) due to the coarse discretization and homogenized material characteristics. Figure 15 shows energy evolution with respect to time for each computational model. The total energy (E t ) consists of strain energy (E s ), kinetic energy (E k ), and fracture energy (E f ). Energy evolution of the adaptive model is similar to the evolution of the direct model, although slight differences are observed.

EFFECT OF MICROSTRUCTURE
Based on the proposed computational framework, the effects of microstructure on the dynamic cohesive fracture behaviors are investigated. The single edge notched tension example in Section 4.3 is employed; in this case, the interfacial bonding strength and unit cell arrangement are varied.

Effect of the interfacial bonding strength
The effect of the interfacial bonding strength on dynamic fracture behavior is investigated with respect to the loading rate. are 600 MPa and 78.5 N/m, respectively. For the weak bonding case, the cohesive strength is 300 MPa, and the fracture energy is 39.2 N/m. Other material properties of the matrix and particles are the same as the material properties described in Section 4.3. Because of the randomness associated with particle position and finite element discretization, 10 computational simulations are performed for each case to evaluate the average crack velocity and dissipated energy, that is, 60 simulations in total. The crack patterns corresponding to the weak and strong bonding cases are illustrated in Figures 16 and 17, respectively. In the weak bonding case, cracks are generally initiated from the interface between the matrix and particles due to the relatively low cohesive strength. Later, some interfacial cracks coalesce while a main crack propagates, and others are arrested in the matrix. In the case of the strong interface, as shown in Figure 17, cracks generally propagate along the matrix and few interfacial cracks are observed. Some of the branching cracks are arrested by the particles. Additionally, when the boundary velocity increases, more local micro cracks are observed for both weak and strong bonding cases. The weak bonding case ( Figure 16C) induces longer microbranching cracks than the strong bonding case ( Figure 17C), especially at a higher boundary velocity, for example, 250 mm/s. The average crack velocities of the weak and strong bonding cases are evaluated according to the boundary velocity, as shown in Figures 18A,B,  tip position of 2.5 mm to a crack tip position of 5.5 mm. Because multiple cracks generally appear during dynamic fracture computations, a crack tip position is defined as the farthest position from the initial notch where a complete failure condition is observed. The average crack speed of the weak bonding case is slightly higher than the speed of the strong bonding case. This is because cracks mainly propagate along the matrix for the strong bonding case, while cracks are generally initiated from the interfaces in the weak bonding case, as observed previously. Thus, early formation of interfacial cracks can increase average crack velocity in the weak bonding case. Figure 19 illustrates the changes in dissipated fracture energy with respect to the changes in the boundary velocity. One notes that a higher boundary velocity corresponds to a higher loading rate. An increase in the loading rate results in an increase in the dissipated fracture energy. This is because a higher loading rate results in more microcracks, as shown in Figures 16 and 17, and thus a greater amount of dissipated fracture energy. Furthermore, such behaviors may be associated with the rate effects of materials such that the strength of a material increases with an increase in the strain rate. 59 This is because more energy is dissipated at a higher loading rate, which can lead to a higher load bearing capacity. According to the increase in loading rate, the increase in dissipated fracture energy for the weak bonding case is more significant than that for strong bonding. This is because more and longer crack branches are observed in weak bonding cases than in strong bonding cases while the boundary velocity increases.

Effect of the unit cell arrangement
To investigate the effect of the unit cell arrangement, two types of unit cell structures are employed, as shown in Figure 20. respectively. An initial boundary velocity of 250 mm/s is applied, and the weak and strong interfacial bonding strengths are used, as defined in Section 5.1.
In the strong bonding case, the crack pattern of unit cell type I ( Figure 21A) is different from the crack pattern of unit cell type II ( Figure 21B). In unit cell type I, diagonal microcracks are arrested by particles and thus stop propagating. In unit cell type II, diagonal microcracks propagate between particles. Therefore, the length of microcrack branching in unit cell type I is generally shorter than that in unit cell type II, while a main crack propagates along the horizontal direction in both types.
However, in the weak bonding case, both unit cell types exhibit similar crack patterns, as shown in Figure 22. Multiple cracks are initiated along the interface between the matrix and particles due to the weak bonding strength, and randomly distributed branching cracks are observed in unit cell types I and II. Although diagonal microcracks are blocked by particles in unit cell type I, they keep propagating along the interface between the matrix and particle due to the weak bonding strength.

CONCLUSION
A computational framework for multiscale dynamic fracture analysis is proposed by integrating far and local fields in conjunction with adaptive microstructure representation. The microstructure of a material is explicitly and adaptively represented at a local field (eg, potential hot spots), while being homogenized using the Mori-Tanaka method at a far field. Then, the far and local fields are solved at the same time using the equation of motion, while physical field quantities such as velocity and displacement are consistently updated by performing local analysis of a unit cell. Additionally, in this study, 4k meshes are simply employed to adaptively refine a local field. The proposed adaptive multiscale analysis framework in not limited to 4k meshes, but can be utilized other meshing techniques for the local field representation. The computational framework is verified by checking energy balances and comparing the proposed model with direct computation. Adaptive microstructure representation results in a kinetic energy difference of 2% between the homogenized material and material with microstructure owing to material-characteristic transitions. A strain energy difference of less than 1% is observed during the local analysis of a unit cell. One notes that the local analysis of a unit cell provides more accurate stress fields than a simple interpolation approach. Additionally, elastic wave propagation and dynamic fracture examples are investigated using adaptive microstructure representation and direct computation. The global displacement, crack patterns, and energy evolution of the adaptive modeling agree well with those of direct computation; furthermore, the computational cost of the proposed approach is significantly lower than that of direct computation. This is because local stress concentrations around crack tip regions are efficiently and effectively captured in adaptive microstructure representation.
Furthermore, the effects of the microstructure of a material on nonlinear dynamic fracture phenomena are investigated. A potential-based cohesive zone model is utilized to approximate the nonlinear fracture process zone, while discontinuities like cracks are explicitly described by adaptively introducing surface elements in conjunction with a topology-based data structure. Subsequently, the effects of interfacial bonding strength, loading rate, and unit cell arrangement are explored on a two-phase ceramic composite material via a single edge-notched tension test. Computational results illustrate that the overall length of microbranching increases with a reduction in bonding strength and/or an increase in the loading rate and thus the overall crack velocity and dissipated energy increase. Additionally, the unit cell arrangement can influence dynamic fracture behavior because microcracks in the matrix can be arrested by particles with a strong bonding strength. In summary, local crack patterns change according to the combination of unit cell arrangement and interfacial bonding strength; therefore, the microstructure of a material should be carefully considered in dynamic cohesive fracture investigation.