Lattice Approach in Continuum and Fracture Mechanics

Lattices have so far been considered in solid mechanics mostly as numerical tools to analyze crack path and geometry in boundary value problems where strong discontinuities are imminent [1–8]. The fundamental concept of using the lattice formulations in solid mechanics was first introduced by Hrennikoff [1] in 1941. He replaced a continuum elastic domain by a lattice of truss, or frame elements, the elastic properties of which were determined based on those of the continuum domain with no cracking. Herrmann [2], in the context of theoretical physics, discussed that a material can be modeled as a lattice of brittle breaking beam elements. Schlangen and van Mier [4] used triangular lattice mesh projected on the grain skeleton to simulate the crack path propagation in heterogeneous concrete materials. They considered three different phases including aggregate, matrix, and their bond and analyzed the branching and bridging with respect to bond strength. Schlangen and Garboczi [5], also, employed regular and random lattices to simulate the geometry and crack pattern of simple shear experiment on a concrete panel. They evaluated three types of elements, i.e., trusses with one and two, and frames with three degrees-of-freedom at each node. Schlangen and Garboczi concluded that in the simulations with frame elements, the crack pattern on the mesh was much closer to the experimentally obtained cracks. Bolander and Saito [6] used an innovative method to discretize the continuum domain, resolving the problem of assigning appropriate values for each strut's cross section area and moment of inertia in the lattice. They used Voronoi polygons to discretize the continuum domain into an assemblage of convex rigid particles interconnected along their boundaries through flexible common sides or interfaces. They established the stiffness matrix for each strut based on the elastic properties of the continuum material and the frame's geometry. van Mier [7], in his book, thoroughly discussed the application of lattice models in concrete fracture including failure aspects of the lattice and its mesh generation. Using the mechanics of materials approach, he compared lattice simulations in terms of crack path and shape with experimental results for concrete specimens under different loading conditions. Tankasala et al. [8] explored crack-tip fields including plastic zone size and crack opening displacement along with the fracture toughness of 2D elastoplastic lattices for four topologies. In their study, it was revealed that the fracture toughness of ductile lattices is sensitive to the length scale of lattice in addition to relative density and choice of topology.

In all the above studies, the main attention was paid on the global behavior of the problem in terms of load–displacement or load–crack mouth opening displacement diagrams along with the crack path and geometry. Mohammadipour and Willam [9] showed that it is possible to evaluate fracture properties of material with evolving cracks like energy release rate, loading phase angle, and stress intensity factors by using the lattice with a conceptually simple linear elastic erosion algorithm. In this study, it is attempted to demonstrate the fundamental capabilities of the lattice approach to simulate the elastic and fracture-based behavior of a 2D plane stress rectangular panel with a single circular hole under different loading conditions, as a classic problem in LEFMs.

Section 2 briefly explains the implementation of the lattice approach, its theoretical framework, geometry, constitutive mechanics, and the governing fracture criteria used for each lattice member's erosion. In Sec. 3, it is shown that how the discrete lattice approach may be used in obtaining the continuum strain distributions and classic SCFs of the aforementioned problem under different loading scenarios. Section 4 deals with extracting the fracture properties of the same problem under direct tension when the crack initiates and emanates from the circular hole.

Hrennikoff [1] in 1941 first introduced how the main concepts of lattices could be used in discretizing a continuum domain and in analyzing its response in terms of applied loads. Voronoi diagrams, based on a nonrandom distribution of points, were used in this study to discretize the continuum domain into an assemblage of convex rigid particles interconnected along their boundaries through flexible common sides or interfaces. Since the example problem here is a rectangular domain with a circular hole, special care is exercised to distribute the computational points over the rectangular domain with smooth transition from polar to Cartesian coordinates. This will be explained in Sec. 3. For more details on the lattice geometry and Voronoi diagram, refer to Refs. [6,9], and [10].

One of the attractions of lattice approximations is the combination of the mechanics model and the material structure, called material structure overlay. The lattice is the mechanical model; the material structure is simply projected on top of the lattice and various properties are assigned to the lattice elements depending on their specific location [6,7,9,10]. In other words, the lattice and material structure are two independent features of the model.

Furthermore, 2D frame elements with three degrees-of-freedom were considered for the numerical lattice simulations in this study. As shown in Fig. 1, each frame strut can transfer, in general, normal force N, shear force Q, and bending moment M, corresponding to the degrees-of-freedom at the nodes.

The relation between these forces and their corresponding displacements at the endpoints of the frame element can be found in the literature [5]. An approach, introduced by Bolander and Saito [6], was employed here to establish the constitutive relation of the continuum in a 2D plane stress lattice simulation. Figure 2 shows the relationship of two Voronoi triangular particles by putting a flexible interface between them. In this figure, two triangular particles are connected at their interfaces by translational and rotational stiffnesses, i.e., $kn$⁠, $kt$⁠, and $kϕ$⁠, which approximate the elastic properties of the continuum. Points 1 and 2 are the nuclei or computational points lying at the beginning and end points of the strut or frame element. Further details are elaborated in Ref. [6].

The simulation of fracture was performed with a “linear elastic” analysis of the lattice under loading, Fig. 1(b), and removing an element from the mesh which exceeds a certain fracture criterion, for instance, a tensile or compressive stress based on the failure envelope. The “gap” between the remaining elements is considered as a discontinuity or crack in the lattice mesh. After removing the element, the lattice mesh contains one less element. The simulation is continued by performing a linear elastic analysis of the new configuration, where the forces that were carried by the removed element are now redistributed over the neighboring elements. This procedure continues until the next element satisfies its “fracture criterion” and so on [5,7]. Thus at each step, the external load on the lattice is increased and the critical element at the fracture threshold is removed. The erosion strategy leads to an “instantaneous relaxation” of the load, carried by that removed member of the lattice [7]. This was observed during the lattice analyses of this study as a sudden drop in form of snapback in the load–displacement diagrams.

where $N$ is the normal force of the lattice element; $A$ is its cross-sectional area; $Mi$ and $Mj$ are the bending moments at the nodes $i$ and $j$⁠, respectively; $Sa$ is the section modulus; $ft$ and $fc$ are the tensile and compressive strength values of the material, respectively; and $0≤α′≤1$ is added to limit the effect of bending in the fracture law. The value of $α′$ may be determined by parametric studies and comparison with experimental measurements [4].

As $n$ increases, the shape of curve approaches to a rectangle. For the case of a circular hole, $â=b̂$ in Eq. (3).

Figure 3 illustrates the lattice mesh generation of the rectangular domain with a circular hole of radius $0.25 (in.)(or 0.64 (cm))$⁠. The domain is a $4(in.)×4(in.)×0.25(in.)(or 10.16(cm)×10.16(cm)×0.64(cm))$ panel, the material properties of which are selected from the brick properties used previously by the authors using a digital image correlation system, which is used in a wide variety of applications from fracture to damage mechanics (⁠$E=4,506,000 (lbf/in.2)(or 31067.78 MPa), ν=0.1$⁠) [12–17]. The tensile of compressive strength of the brittle brick is $400(lbf/in.2)(or2.76MPa)$ and $6000(lbf/in.2)(or41.37MPa)$⁠, respectively. Due to the symmetry along $x$ and $y$ axes, one-quarter of the panel is analyzed, for which the mesh has number of degrees-of-freedom and elements of 2776 and 2594, respectively. Figure 3(a) shows the distribution of centroids, or computational points, over the panel using Eq. (3) with $â=b̂$⁠. Voronoi particles are then defined based on these centroids. In Fig. 3(b), lattice struts or frame elements were established through connecting the centroids of neighboring particles in Fig. 3(a). The parameters used to define the lattice mesh are $S$⁠, i.e., $Smin$⁠, $Smax$⁠, and $Δθ$⁠. $Smin$ and $Smax$ are the smallest and largest values of $S$ creating an adaptive mesh which take into consideration the strain gradients of the domain, which are the highest close to the hole. Depending on the size of the hole, $R$⁠, and the width of the panel, $2b$⁠, the parameter $S$ should be determined from Fig. 4 for a selected value of $Δθ=π/64$⁠. This curve was obtained by the lattice in uniaxial tension where a sensitivity analysis was conducted for a range of $R$ and $b$ values. The objective was to obtain an SCF of 3 with a maximum error of $5%$ while changing the above mentioned mesh parameters.

A linear elastic analysis was conducted to obtain the strain distributions due to the relative displacements of each computational point at the lattice. The strain distribution calculation out of the lattice analysis is a postprocess stage at the end of each increment during the linear elastic behavior of the panel. The strain distribution plots are like an extra layer on top of the lattice mesh. To upscale the mesolevel lattice results into continuum strain distributions, an imaginary continuum Q4 finite element is mapped on top of the lattice mesh, the nodes of which exactly coincide with neighboring computational points of the lattice mesh (Fig. 5).

It is seen in Eq. (5) that the effect of rotation is only present at shear strain components of $ε$⁠, and the normal strains are similar to those of the classic continuum mechanics. Similar to the classic argument of calculating gradients at Gauss points in the finite element formulation, the components of $ε$ are obtained at one Gauss point or centroid of each Q4 element, given the translational and rotational nodal values obtained from the lattice.

An extrapolation technique is needed to determine the strain values at nodes of each Q4 element to have those values at the interior side of the circular hole. Zienkiewicz et al. [21] developed a simple method named superconvergent patch recovery (SPR) for determining the derivatives (strains and stresses) of the finite element solutions at nodes. In this method, the derivative values at Gauss points inside an element are smoothed by a polynomial of order $p$ within a patch of elements for which the number of sampling point (Gauss points) is greater than the number of parameters in the polynomial. The strain in the patch is expressed in global coordinates by a polynomial expression, the parameters of which are determined by a least square method. The SPR method was used in this study to extrapolate the strain values at nodes. It should be noted that the gradient values obtained by the SPR method at internal nodes of a regular mesh are the same as those obtained by averaging over the neighboring elements. However, this is not true when the mesh is irregular or boundary nodes are of interest like nodes $A$ and $B$ in Fig. 6(a). Therefore, a proper extrapolation technique like SPR will retain high accuracy [21].

The lattice mesh with circular hole in Fig. 3(b) is considered under three loading scenarios, namely, uniaxial tension in $y$ direction, biaxial tension in $x$ and $y$ directions, and biaxial tension–compression in $y$ and $x$ directions, respectively. For an infinite thin element in uniaxial tension with a single circular hole, the SCF for $εy$ is theoretically $3$ at point A shown in Fig. 6(a). It should be noted that the value of $3$ belongs to classic linear elasticity where there is no rotational degree-of-freedom as opposed to the micropolar elasticity that the rotations have reduction effects on the concentration factors around the hole [22]. The SCF for $εy$ and $εx$ at points A and B for the same panel with biaxial tension in $x$ and $y$ directions is $2$⁠, while for the case of biaxial tension–compression in $y$ and $x$ directions, the SCFs at points A and B are $4$ and $−4$ for $εy$ and $εx$⁠, respectively [23]. Figure 6(b) illustrates the distribution of SCF for $εy$ for the panel in uniaxial tension. The SCF evaluated by the lattice approximation at point A is $2.97$⁠. Figures 6(c), 6(d), 7(c), and 7(d) show the distribution of SCF for $εy$ and $εx$ for the biaxial tension–tension and tension-compression, respectively. In Figs. 6(c) and 7(c), the SCF for $εy$ and $εx$ at points A and B is equal to $1.88$⁠, respectively. For the case of biaxial tension–compression which is pure shear, SCFs assessed by the lattice are $3.94$ and $−3.94$ at points A and B for $εy$ and $εx$⁠, respectively. Table 1 compares the SCF values of the lattice and those obtained from the classical theory of elasticity for the problem in Figs. 6 and 7. It is seen that the discrete microlevel lattice formulation can well analyze the SCF of a 2D plane stress panel in the aforementioned loading conditions. The lower obtained values by the lattice compared to those of the classic theory of elasticity are due to the effect of the rotational degrees-of-freedom in the lattice which increase the stiffness of the mesh, thus resulting in lower SCFs. In this section, it was shown that the lattice results may be considered in the continuum sense to determine strain and stress distributions of a domain in the context of linear elasticity.

The implemented 2D lattice explained in Sec. 2 is capable of simulating crack path evolution in the form of strong discontinuities at a homogeneous or heterogeneous solid. Since the crack propagation is captured by the lattice during an analysis, it is postulated that the fracture mechanics quantities like the energy release rate or the stress intensity factors associated with the evolving crack may be determined by the lattice. Other numerical techniques like the classic finite element method of virtual crack extension [24–27] or the well-known extended finite element method (XFEM) [28–30] may also be used for obtaining the fracture quantities of crack problems. In the virtual crack extension procedure, as used by Charalambides et al. [25] and Matos et al. [26], a precracked finite element mesh with length $a$ was considered and the virtual crack extension method was applied by virtually increasing the crack length and changing the stiffness of a ring of elements around the crack tip. This method is not based on a progressive crack evolution where the crack length “$a$” increases during a single simulation. XFEM, which is based on the mathematical foundation of the partition of unity FEM, is also capable of measuring the fracture quantities. However, its implementation needs much more effort and considerations than the lattice regarding, for instance, modeling the arbitrary crack propagation paths handled by “level set” method, multiple crack configurations, cracks intersecting with other discontinuities, and also cracks emanating from holes or other internal interfaces. Although the goal here is not to compare the capabilities of XFEM and the lattice, its relative simplicity both in theory and implementation along with accepted results [9] make the lattice approach an interesting approach in solving crack problems in the context of fracture mechanics.

where $[K]a+Δa$ is the stiffness matrix after the crack growth $Δa$⁠. Equation (6) can be employed by the lattice to determine the energy release rate of a propagating crack, merely by an energy method and change in the configuration of the lattice mesh and its computational points.

where $Δu$ and $Δv$ are the relative crack flank displacements of two points on the top and bottom of crack surfaces at a distance $r$ behind the crack tip defined as $Δu=u(r,θ=π)−u(r,θ=−π)$⁠, $Δv=v(r,θ=π)−v(r,θ=−π)$ (Fig. 8). Values of $Δu$ and $Δv$ are directly obtained from the lattice simulation as the crack propagates at a distance $r$ from crack tip to the location where the most recent strut has just been removed upon failure. The stress intensity factors for modes I and II are then determined using the phase angle by

Newman in 1971 [32] analytically solved the problem with a crack emanating from a circular hole in a rectangular panel subjected to uniaxial tension, as shown in Fig. 9 [33]. He expressed the mode I stress intensity factor as

where $σ$ is the far-field tensile stress; $a, b, R, and h$ are illustrated in Fig. 9; and $F$ is a function of crack length, panel, and perforation geometry determined analytically to modify the value of $K1$⁠. Newman obtained $F$ for a panel with $h/b=2$ and $R/b=0.25 and 0.5$⁠.

A 2D plane stress rectangular panel with the dimensions of $8(in.)×4(in.)×0.25(in.)(or 20.32(cm)×10.16(cm)×0.64(cm))$⁠, with a circular perforation of radius $1 (in.)(or 2.54 (cm))$⁠, and with the same material properties mentioned in Sec. 3 was analyzed by the lattice formulation in uniaxial tension. Because of symmetry, a quarter of the panel was considered. Thus, for this analysis, $h/b=2$ and $R/b=0.5.$ Figure 10(a) illustrates the lattice mesh used to simulate the Newman problem shown in Fig. 9. Figure 10(b) belongs to the middle of lattice analysis when the crack propagated form the circular perforation emulating the Newman problem when $a>R$⁠. Figure 11 depicts the lattice simulation results for this problem. Figure 11(a) shows the load–displacement curve of the panel with snap-back instabilities with instantaneous relaxation of the load, carried by that removed part of the lattice. The saw-tooth pattern observed at the postpeak part of the curve is caused by failure of each strut in an unzipping manner, releasing fracture surface energy in the form of displacement relaxations of the lattice mesh. In Fig. 11(b), the analysis results of mode I stress intensity factor, $K1$⁠, have been compared with those of Newman analytical solution in the terms of function $F$ in Eq. (13). The energy release rate was determined by the lattice according to Eq. (6) with considering both displacement vectors ${u}$ before and after crack length increment of $Δa$⁠, i.e., ${u}a$ and ${u}a+Δa$⁠, at each time step as

where $Ga$ and $Ga+Δa$ are obtained using ${u}a$ and ${u}a+Δa$ using Eq. (6), respectively. The energy release rate obtained by the lattice is used as an input in Eq. (9) to obtain the modulus of stress intensity factor. The lattice also provides the crack flank displacements $Δu$ and $Δv$ used in Eq. (10) to determine the loading phase angle which is an essential gradient to obtain $K1$ employing Eq. (11). The loading phase angle has the minimum, maximum, and average values of $0.24deg$⁠, $5deg$⁠, and $2.73deg$⁠, respectively, for $0.51<a/b<0.91$⁠, indicating that mode I is the prominent failure mechanism. $K1$ is then divided by $σπa$ with updated $σ$ and $a$ at each time step to obtain $F$⁠. As it can be seen in Fig. 11(b), the lattice approach can accurately assess the stress intensity factor of the Newman problem with an error which gradually increases from $0.11%$ when $a/b=0.51$ to $10.66%$ when $a/b=0.9$ with an average value of $5.56%$⁠. A possible explanation for this slight increase might be that the lattice adaptive mesh gets coarser in the radial direction of the polar coordinates as the crack evolves and as the crack tip distance from the center of the hole increases (Fig. 10(b)).

Lattice model was investigated from two perspectives in the context of both continuum and fracture mechanics. A well-known classic problem in the theory of elasticity was simulated by the implemented 2D lattice. In the context of Cosserat continuum mechanics and linear elasticity, a rectangular plane stress thin panel with a single circular hole was solved by the lattice approach subjected to uniaxial, biaxial, and shear loadings. The translational and rotational displacements of the lattice's computational points were used as inputs in a Q4 continuum finite element to extract the strain distributions and SCFs around the hole. It was observed that the lattice well evaluates the SCFs for $εy$ and $εx$ in the mentioned loading conditions.

In the context of LEFM, the lattice was employed to derive the fracture quantities of the same problem with a perforation using an energy method. Energy release rate, loading phase angle, and then stress intensity factor for mode I were evaluated. The lattice results of fracture quantities were then validated by comparing with the analytical solutions of Newman for the same problem. Though simple in constitutive formulation and failure criterion, the lattice approach may be used in obtaining these fracture quantities of other problems [9] where there is a discontinuity in the form of a propagating crack.

The authors would like to acknowledge the partial support of this research effort by the U.S. National Science Foundation under the NSF Project No. 1100971 named “Interface Mechanics of Masonry Panels under Bi-axial Loading.” The authors would also like to appreciate the valuable comments of the anonymous reviewer.

$a$ =

crack length

$A$ =

cross-sectional area for a lattice element

$â,b̂$ =

half of an ellipse's minor and major axes, respectively

$b$ =

half of the panel width

$E$ =

Young's modulus

$E*$ =

equivalent Young's moduli for plane stress and plane strain

$F$ =

modification coefficient for $K1$

$ft,fc$ =

tensile and compressive strengths of material

$G$ =

energy release rate

$h$ =

half of the panel's height

$h0$ =

lattice element length

$[K]$ =

global stiffness matrix of the whole lattice mesh

$|K|$ =

modulus of complex stress intensity factor

$kn,kt,kϕ$ =

normal, tangential, and rotational stiffnesses connecting two particles in the lattice formulation

$[K]a+Δa$ =

global stiffness matrix of the whole lattice mesh after a crack growth of $Δa$

$K1,K2$ =

stress intensity factor for modes 1 and 2, respectively

$l43$ =

lattice element width

$M$ =

internal bending moment of the lattice frame element

$Mi,Mj$ =

bending moments of the lattice frame element at the nodes $i$ and $j$

$N$ =

internal normal force of the lattice frame element

$Q$ =

internal shear force of the lattice frame element

$r$ =

radial distance from the crack tip

$R$ =

radius of the circular perforation

$S,Smin,Smax$ =

mesh parameters

$Sa$ =

section modulus

SCF =

strain concentration factor

${u}$ =

displacement vector of the lattice mesh

$u,v$ =

components of displacement along $x$ and $y$ axes, respectively

$x,y$ =

global coordinates

$α′$ =

scalar factor to limit the contribution of bending in the lattice frame's fracture criterion

$α*$ =

bi-elastic coefficient in the Hilbert equation

$δn,δt,ϕ$ =

local displacements in the normal, tangential, and rotational directions

$Δu,Δv$ =

relative crack flank displacements along $x$ and $y$ axes, respectively

$Δθ$ =

mesh parameter

$εij$ =

strain tensor

$εijk$ =

Levi-Civita tensor

$ν$ =

Poisson's ratio

$Π$ =

total potential energy

$σ$ =

normal stress

$σeff$ =

effective stress defined in the lattice frame's fracture criterion

$ψ$ =

loading phase angle