Access the full text.
Sign up today, get DeepDyve free for 14 days.
iitm.sudhakar@gmail.com Institute for Computational We devise a ﬁnite element methodology to trace quasi-static through-thickness crack Mechanics, Technical University paths in nonlinear elastic solids. The main feature of the proposed method is that it can of Munich, Boltzmannstr. 15, 85747 Garching, Germany be directly implemented into existing large scale ﬁnite element solvers with minimal eﬀort. The mesh topology modiﬁcations that are essential in propagating a crack through the ﬁnite element mesh are accomplished by utilizing a combination of a mesh reﬁtting procedure and a nodal releasing approach. The mesh reﬁtting procedure consists of two steps: in the ﬁrst step, the nodes are moved by solving the elastostatic equations without touching the connectivity between the elements; in the next step, if necessary, quadrilateral elements attached to crack tip nodes are split into triangular elements. This splitting of elements allows the straightforward modiﬁcation of element connectivity locally, and is a key step to preserve the quality of the mesh throughout the simulation. All the geometry related operations required for crack propagation are addressed in detail with full emphasis on computer implementation. Solving several examples involving single and multiple cracks, and comparing them with experimental or other numerical approaches indicate that the proposed method captures crack paths accurately. Keywords: Mesh reﬁtting method, Crack propagation in nonlinear materials, J-integral, Nodal releasing technique, Fracture of nonlinear solids Background Objective and motivation One of the main reasons why devising a computational methodology to deal with fracture mechanics is challenging, is the fact that cracks propagate in arbitrary directions through the material. If the dynamics of the crack were known apriori, one can design an optimal mesh that allows the propagation of a crack through the pre-existing mesh at each instant. Since this is not the usual case, the mesh has to be repeatedly modiﬁed to accommodate the advancement of cracks within the ﬁnite element (FE) mesh. The objective of this work is to devise a simple procedure to achieve the required mesh modiﬁcations, which enables us to model complex crack paths through nonlinear elastic solids. The present work is motivated by our interest in developing computational methodologies for ﬂuid- structure-fracture interaction (FSFI) [1] that model the following phenomenon: when a ﬂexible structure interacts with the ﬂuid ﬂow, the ﬂuid loading induces elastic deformation © The Author(s) 2017. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. 0123456789().,–: vol Sudhakar and Wall Adv. Model. and Simul. in Eng. Sci. (2017) 4:2 Page 2 of 23 as well as fracture failure of the structure, and the ﬂuid medium ﬁlls the crack opening. The ﬁrst step of extending ﬂuid-structure-interaction (FSI) methods to handle FSFI is to equip the structural analysis with a fracture mechanics solver. This is achieved in this work by developing a crack propagation approach which facilitates, with minimal implementation eﬀorts, • to update the existing large scale structural mechanics solver into a robust tool to handle single and multiple quasi-static cracks • to couple the crack propagation method with existing FSI approach to model FSFI. The method devised in this paper is implemented in BACI, a large scale parallel multi- physics solver developed at our institute. It is to be mentioned that the objective of the present work is not to devise a method which is competitive to the available class of meth- ods in terms of computational eﬃciency or accuracy. Rather, the focus is on devising a simple crack propagation approach, which circumvents the complexities associated with them (as discussed below) to aid the development of an FSFI solver. Brief overview of relevant methods The present work is based on fracture mechanics [2–12] framework rather than on contin- uum damage mechanics [13,14] because we model the propagation of sharp cracks in the material. Majority of the existing computational methods employing fracture mechanics principles utilize either one of the following frameworks: • Adaptive remeshing • Enriched partition of unity These methods are very successful and ﬁnd plethora of applications but implementing them in an existing (large scale) FE package pose several challenges. Adaptive remeshing methods These methods, as the name implies, adaptively reﬁne the mesh in the crack tip vicinity where the solution dictates the dynamics of crack propagation, and coarsen the mesh away from the crack tip. They make use of special data structures [2,3,15], together with either a globally adaptive remeshing procedure [4–6] or a local mesh modiﬁcation algorithm [7]to accommodate crack propagation at arbitrary directions within the computational domain. Each time when a crack extends, these methods introduce several new nodes into the mesh. As a result, they require mesh generation related algorithms to modify the mesh appropriately. Usually in a large scale FE code which is generalized to address multiscale and multiphysics problems, such fracture-speciﬁc and mesh-modiﬁcation routines are neither available nor easy to implement in a generalized way. Enriched partition of unity methods (EPUM) These methods represent the recent developments in computational fracture mechanics. They include extended ﬁnite element methods (XFEM) [16–19] and generalized ﬁnite ele- ment methods (GFEM) [20,21], and these class of methods are originally developed with an objective to eliminate the adaptive remeshing and its associated complex and time- consuming operations. The fundamental idea behind this method is to enrich the ﬁnite element solution space with additional problem-speciﬁc enrichment functions. The crack Sudhakar and Wall Adv. Model. and Simul. in Eng. Sci. (2017) 4:2 Page 3 of 23 can propagate within the interior of an element, and hence it is possible to simulate crack propagation without modifying the underlying discretization. Though these methods are demonstrated to be powerful, the following points hamper their easier implementation into an existing structural mechanics solver. The numerical integration of singular enrich- ment functions is still an active area of research in EPUM [22–28]. Moreover, these meth- ods require the implementation of complicated geometry-mesh intersections for which robustness is always an issue, especially in 3D. Also, the number of degrees of freedom attached to a few nodes changes each time when the crack advances. Most importantly, except a few (e.g. [29–31]), all studies are focused on crack propagation through linear elastic materials mainly because the formulation of enrichment functions in nonlinear regime is still an active area of research. Owing to the aforementioned implementation issues, neither EPUM nor adaptive remeshing methods are ideal for developing ﬂuid-structure-fracture interaction meth- ods. This is because FSFI methods have to handle the combination of challenges from two sources: those arising from crack propagation besides another big challenge from FSI. The aim of this work is to devise a simple crack propagation approach that is suitable speciﬁcally for developing FSFI. The method developed in this work, as will be explained later, shares a similarity with arbitrary Lagrangian Eulerian (ALE) based methods that it involves a mesh-deformation step. Therefore, ALE based methods for fracture mechanics are brieﬂy recalled. The use of ALE in computational fracture mechanics is not widespread. Only few studies employed ALE to address crack propagation problems, and a brief account of majority of such works is provided below. Existing ALE based crack propagation methods Though ALE formulations are widely used in several solid mechanics applications (refer to [32] for an overview), less than a handful of researchers used them to handle fracture mechanics problems. The ﬁrst use of ALE is described in [33,34] to model dynamic crack propagation. The capability of these methods are demonstrated by simulating a few mode- I dynamic crack propagation problems, and comparing them with analytical relations. In this work, the material separation is not explicitly modeled. Moreover an existing Lagrangian FE code framework cannot be directly extended to include this model. In order to achieve this, an improved method has been developed in [35], which is applied to simulate mixed-mode dynamic crack propagation examples. Another method that uses a mesh motion algorithm based on an isoparametric mapping is presented in [36], in which a self-similar dynamic crack propagation problem in double cantilever beam is solved. They conclude that ALE methods are more robust, and can be a powerful alternative to remeshing. For a better understanding of dynamic crack growth in Fibre-reinforced plastic (FRP) composites, an ALE based method together with a contact mechanics approach is developed in [37]. A very ﬁne mesh in the neighborhood of the crack tip is maintained throughout the simulation with the help of a remeshing procedure. This method is then used to study the interfacial debonding phenomenon in FRP strengthened reinforced concrete beams. Another attractive method that combines the advantages of element free Galerkin (EFG) method and ALE is presented in [38]. This method moves a cloud of nodes along with the crack tip, so that the vicinity of the crack tip is always adequately resolved. As EFG is a meshless method that does not require nodal connectivities, such implementation Sudhakar and Wall Adv. Model. and Simul. in Eng. Sci. (2017) 4:2 Page 4 of 23 of ALE to maintain high nodal density in the preferred region is accomplished eﬀectively without resorting to remeshing strategies (refer to [39] for complications involved in implementing such a method in mesh-based FEM). The method is successfully applied to simulate wave propagation and dynamic crack propagation. To summarize, to our knowledge neither a complex trajectory of a single crack nor simple propagation of multiple cracks within a material is modeled until now using ALE based methods. This is predominantly due to the fact that the mesh modiﬁcation method used in ALE, in its classical sense, cannot handle the mesh topology changes that are introduced by advancing a crack through the FE mesh. Continuous remeshing is manda- tory to eliminate the associated mesh tangling problems. In this work, this is avoided by using an additional step in the mesh reﬁtting procedure that allows to modify the element connectivity locally to preserve the quality of mesh, as will be explained in the later section. Structure of the paper The remainder of the paper is structured as follows. The next section briefs the govern- ing equations and the boundary conditions for the problem. Then the complete crack propagation algorithm together with all the details necessary for computer implementa- tion are presented. Finally several numerical examples to demonstrate the accuracy of the proposed method are described. Governing equations s s At reference time t = t , let the structure occupy the domain ,with denoting the boundary of the structure (Fig. 1). is divided into three non-overlapping portions such s s s s s s that = ∪ ∪ in which and are Dirichlet and Neumann portions of D N c D N the boundary respectively, and denotes the crack surfaces which contain always two s s s physical crack faces = ∪ . The balance of linear momentum equation is written c c+ c− as, s s s s ρ d − Div (FS) = ρ b in × (0,T)(1) s c Fig. 1 A schematic representation of a structural domain containing a crack Sudhakar and Wall Adv. Model. and Simul. in Eng. Sci. (2017) 4:2 Page 5 of 23 where ρ is the density of the structure, F is the deformation gradient, S is the second Piola-Kirchhoﬀ stress tensor, b represents externally applied body force per unit mass, 2 s s d d d is the structural displacement and d = . T is the end time of the considered dt time interval, and Div(·) is the divergence operator deﬁned with respect to the material reference frame. Since this is an evolutionary problem involving second order time derivative, initial s dd conditions must be speciﬁed on d and its ﬁrst derivative d = dt s s s s s s ˙ ˙ d | = d on ; d | = d on (2) t=0 t=0 0 0 0 0 Over the boundary of the domain, Dirichlet conditions are speciﬁed on , Neumann conditions are prescribed on , and the crack surfaces are assumed to be traction-free. s s s s s s ¯ ¯ d = d on × (0,T); FS · n = h on × (0,T)(3a) ( ) D N c+ s c− s (FS) · n =0on × (0,T); (FS) · n =0on × (0,T)(3b) c+ c− We deal with hyperelastic Neo-Hookean materials in this work. The strain energy function for such material is given as s s μ λ s 2 = (tr C − 3) − μ ln J + (ln J) (4) NH 2 2 where C is the Cauchy–Green tensor, J is the determinant of deformation gradient J = s s det F, λ and μ are Lame’s constants. The mesh reﬁtting approach The complete numerical methodology, together with the computer implementation aspects, of the present approach are presented in this section. It is assumed that the fracture behavior of the material is completely characterized by the J-integral. The focus of the present work is to simulate through-thickness mixed-mode quasi-static crack propagation within a structure. The current work can be considered as an extension of Tabiei and Wu [40] which describes the implementation of a crack module in DYNA3D FE package, and shares similarities with the method of Miehe and Gürses [3]. Both works address crack propagation through linear elastic materials only. Moreover the complex geometry related operations, like deciding the new crack tip nodes, are not addressed in depth. These details are crucial for implementation of the method. The approach proposed in [3] was called an r-adaptive method, but in order to avoid confusion with complex r- adaptive mesh redistribution methods [41,42], the present method is labelled as mesh reﬁtting approach. The following section presents the complete implementation details of the present method. Though we simulate through-thickness cracks, for simplicity we explain the geometry-related operations in 2D. At each time step, the governing equations are solved by freezing the location of the crack. Then, the solution obtained is used to perform crack propagation related operations as described in Algorithm 1. Sudhakar and Wall Adv. Model. and Simul. in Eng. Sci. (2017) 4:2 Page 6 of 23 Algorithm 1 Computational crack propagation procedure at each time step 1: Construct local coordinate system at crack tip 2: Compute vector J-integral 3: Check whether crack propagates or not, from crack propagation criterion 4: if Crack propagation criterion is not satisﬁed then 5: continue to the next time step 6: end if 7: Obtain the direction of crack propagation from crack kinking criterion 8: Find the new crack tip nodes 9: Apply mesh reﬁtting procedure 10: Propagate the crack using nodal releasing technique Solve the governing equations The ﬁrst step is to solve the structural dynamic equations by freezing the location of the crack. The strong form given in Eq. (1) is multiplied by appropriate test functions (δd ) and are integrated over the structural domain to obtain the weak form which is stated as, s s Find d ∈ W such that for all δd ∈ V , the following holds d d s s s s s s s s s ¨ ¯ s s s (δd , ρ d ) + (Grad δd , FS) s = (δd , ρ b ) +δd , h (5) 0 0 N s s where (., .) and ., . mean the standard L -inner product over the reference domain 0 N and Neumann part of the boundary, respectively. The solution space and the test function space are deﬁned as s 1 s s s W ={d ∈ H ( ) | d = d on } (6) 0 D s 1 s s s V ={δd ∈ H ( ) | δd = 0on } (7) 0 D The above integral equations are dealt with nonlinear FEM for spatial discretization and Generalized-α method for time discretization. The resulting nonlinear system of algebraic equations are solved using Newton–Raphson method to obtain the solution of displacement ﬁeld (d ). For a more elaborate discussion of this well established procedure, the reader can refer to the literature (e.g. [43,44]). Perform computational crack propagation procedure The displacement solution obtained from the previous step is used to perform crack propagation procedure by computing vector J-integral. It involves seven discrete steps, each of which are detailed below. Step 1: Construct local coordinate system at crack tip To compute fracture mechanics quantities from the FE solution, and to decompose these quantities into their corresponding modes in a mixed-mode problem, it is essential to construct a local coordinate system (ξ, η)atthe cracktip (x ) as shown in Fig. 2.The base vectors (e , e ) associated with (ξ, η) can be easily constructed because e is the symmetry 1 2 1 line of the crack; e can be obtained by computing the normal to e in a right hand 2 1 coordinate system. Step 2: Compute vector J-integral The J-integral quantiﬁes the strength of singularity at the crack tip in nonlinear elastic materials [45]. Moreover, it is a single parameter that dictates whether the crack propa- Sudhakar and Wall Adv. Model. and Simul. in Eng. Sci. (2017) 4:2 Page 7 of 23 η, e ξ, e Fig. 2 Construction of local coordinate system at crack tip. The shaded Quads represent ﬁnite elements gates or not, and if at all it propagates in which direction it advances. Hence, it is essential to accurately evaluate the J-integral. With respect to the spatial conﬁguration, the energy release rate along the direction of a crack is deﬁned as, ∂d J = wn − n · σ · dγ (8) ∂ξ where γ is the integration contour, w is the strain energy stored per unit deformed volume, n is the normal to contour γ and σ is the Cauchy stress tensor. All these quantities are deﬁned in the current spatial conﬁguration as shown in Fig. 3a. Since the current work ab Fig. 3 J-integral: a notation, b distribution of support function direction of propagation Sudhakar and Wall Adv. Model. and Simul. in Eng. Sci. (2017) 4:2 Page 8 of 23 employs a total Lagrangian formulation, it is convenient to express the J-integral in the reference conﬁguration using pull-back operations. As shown in [30], it reads as ∂d J = WN − N · P · d (9) where ,W, N are the corresponding quantities in reference conﬁguration, P is the ﬁrst Piola-Kirchhoﬀ stress tensor and ( ,H) denote the crack tip coordinate system in the reference conﬁguration, which is given by its base vectors (E , E )inFig. 3a. 1 2 In FEM, the contour integrals are cumbersome to implement because the material variables are available only at the Gauss points. Interpolating these variables over the desired contour presents complications, in addition to introducing interpolation errors. Hence, several studies [46,47] have proposed the idea of converting the contour integral into an integral evaluated over a ﬁnite domain around the crack tip by applying the divergence theorem. This procedure is straightforward to implement, as it requires only the quantities at Gauss points of elements within the ﬁnite domain. The equivalent domain form of J-integral for large deformation problems is given as [29,30], ∂d J = · P − W I · ∇ (q) dS (10) ∂X where S is the domain enclosed by ,and ∇ (q) is the gradient of support function q with j 0 respect to the reference conﬁguration. In order to construct q, using nodal connectivity information, all the elements that are located on n−layers around the crack tip (see Fig. 3bwith n = 4) are considered. From this, the elements connected to the crack tip are deleted. Then, all the nodes that are on the outer boundary of this element set are located, and among these nodes, the one which has the shortest distance (r ) from the crack tip is chosen. Then the support function is min initialized to take a value of unity at the inner layer of nodes, and drops smoothly to zero when the distance of a node from crack tip is more than or equal to r . The distribution min of q within the integration domain is given in Fig. 3b. Note In this work, we assume that the crack surfaces are traction-free. However, in FSFI applications, ﬂuid loads are acting on the crack faces, and as a result an additional term appear in the computation of J-integral. This is explained further in [1]. Step 3: Check crack propagation criterion The crack propagation criterion determines whether the existing crack propagates through the structure under the current stress state. Crack propagation occurs when the driving force reaches or exceeds the material resistance. The J-integral provides a measure of driving force, and its critical value (J ) is assumed to be a material property, which quantiﬁes the material’s resistance to crack propagation. The present work makes use of the vector J-integral based crack propagation criterion proposed by Ma and Korsunsky [48]. This criterion requires, ﬁrst, the calculation of maximum strain energy release rate, which is given by the magnitude of J. 2 2 G = J + J (11) max 1 2 Sudhakar and Wall Adv. Model. and Simul. in Eng. Sci. (2017) 4:2 Page 9 of 23 where J and J are the strain energy release rates along the direction of e and e respec- 1 2 1 2 tively, which are given by simple dot products J = J · e and J = J · e . 1 1 2 2 The crack extension under the given loading conditions occurs, when G reaches the max fracture toughness of the material (J ), G ≥ J (12) max c If the crack propagation criterion is not satisﬁed, then there is no need to perform the remaining operations; the algorithm moves to the next time step to solve the governing equations of the structure. Step 4: Obtain the direction of crack propagation After conﬁrming that the crack propagates, the next logical step is to determine along which direction it is going to advance through the material, which is provided by the crack kinking criterion. There are several methods put forward to determine the crack kinking direction, and the most important methods are • Maximum circumferential stress criterion [49], • Minimum strain energy density criterion [50], • Maximum energy release rate criterion (MERR)[48,51]. It is concluded in a comparative study [52] that the minimum strain energy density criterion is less accurate, and the accuracy of maximum circumferential stress criterion and maximum energy release rate criterion are equivalent in all the tests considered. This work incorporates the maximum energy release rate criterion, proposed in [48], which is consistent with the crack propagation criterion given in the last section. It is stated that the crack propagates when G reaches or exceeds the characteristic fracture max toughness of the material. MERR predicts the crack propagation direction (θ )tobethe direction of J, which is simply given as −1 θ = tan (13) It is to be remembered that θ is measured with respect to the crack normal, as indicated in Fig. 2. Moreover, in this work, the extent of crack propagation is always set to be the length of one complete edge of an element. Step 5: Find new crack tip nodes Having computed the crack propagation direction from J, the next essential step is to determine the new crack tip nodes i.e, nodes in the FE mesh through which the crack must be propagated. A geometry-based method is used in this work to identify the new tip nodes. The ﬁrst step is to identify all the elements that are connected to the current tip node (Fig. 4). Among the edges of these elements, the edge that is intersected by the propagation vector is found. This intersecting edge is drawn with a thick continuous line in Fig. 4. Then the angles formed by the line joining the current tip node to the edge nodes, and the crack propagation direction are calculated (φ and φ in ﬁgure). 1 2 Sudhakar and Wall Adv. Model. and Simul. in Eng. Sci. (2017) 4:2 Page 10 of 23 current tip diagonal nodes propagation vector Fig. 4 Procedure to ﬁnd new crack tip nodes After getting the required intersecting edge, the next step is to check whether the crack propagates along the diagonal. This is realized by the condition, φ ≤diag-tol (=0.25 radians in all the simulations). In this case, the diagonal node corresponding to φ is marked as new tip node. Since the crack propagates through a diagonal of the element, this element must be split along this diagonal to accommodate crack propagation through the mesh, as explained in the next step. If φ >diag-tol, then the non-diagonal node of the intersecting edge will be the next new tip node. In either cases, all the new tip nodes (if there are multiple cracks, each crack tip will have its own new tip node) are stored in R . Moreover, the distance between the ale intersection point and the new tip nodes, marked as δ in the ﬁgure, are computed and ale will be used in the next step. Step 6: Mesh reﬁtting procedure In this step, we reﬁt the existing mesh in such a way that the modiﬁed mesh contains an edge along which the crack can propagate. The mesh reﬁtting procedure used in the present work involves two discrete operations listed as follows: 1. Nodal repositioning 2. Splitting quadrilateral (Quad) elements into triangular (Tri) elements In the ﬁrst step, the nodes are repositioned without touching the elements. This means that the element topology (shape and total number of elements, or the connectivity between elements) remains unchanged. The elements only deform due to the movement of the nodes. Nodal repositioning in this work is achieved by solving elastostatic equations but obvi- ously any other mesh moving approach could be used as well. In this approach, the mesh is treated as a linear elastic body, and the governing equations of the mesh movement together with the boundary conditions are ale 1 Sudhakar and Wall Adv. Model. and Simul. in Eng. Sci. (2017) 4:2 Page 11 of 23 m s ∇ · σ = 0on (14a) d = 0on ∂ (14b) ale d = δ on R (14c) ale ale m m where σ is the ﬁctitious Cauchy stress tensor, d denotes the displacements at each node within the mesh, represents the whole structural domain, and ∂ denotes the ale s s s boundary for ALE computations: ∂ = ∪ ∪ . Displacements at the new crack ale D N c tip nodes are set to be δ that is computed in the previous step. ale The above equations are solved to obtain the mesh displacement d , which is used to move each node in the mesh to its new location. After this mesh movement operation, the new crack tip node is moved to the intersection point along the intersection edge (see Fig. 4). In the next step of the mesh reﬁtting procedure, the Quad elements that are marked to be split are cut into two Tri elements. This happens when the crack propa- gates very close to the diagonal of a Quad element (Fig. 5a).Thisprocessdoesnot involve introducing new nodes into the mesh. By comparing Figs. 4 and 5a, the eﬀect of the mesh reﬁtting procedure is clear: the new tip nodes are ﬁrst moved to the desired location using the nodal repositioning step, and then the Quad element is appropriately split into Tri elements to enable crack propagation along the diago- nal. In short, the combination of nodal repositioning and element splitting ensure that after the mesh modiﬁcations, the crack propagates along an existing edge in the new mesh. One of the main reasons for the failure of ALE based methods in handling large defor- mation or topology change is that such methods maintain their nodal connectivity during the entire simulation. The element splitting operations used in the present work allevi- ates this problem by enabling us to modify the connectivity between the elements locally. This is an essential step without which the nodal repositioning method cannot handle the change in mesh topology that is inherent to crack propagation problems, without resorting to complicated and time-consuming remeshing procedures. Fig. 5 Nodal release technique a splitting an element to allow crack propagation along diagonal, b modify connectivity locally near crack tip; the crack opening is shown only for visualization Sudhakar and Wall Adv. Model. and Simul. in Eng. Sci. (2017) 4:2 Page 12 of 23 Step 7: Nodal releasing technique The two previous steps have enabled us to identify new tip nodes, and to move these nodes to match the computed propagation angle. However, the material separation is not yet included within the FE procedure. In order to achieve this, and to form physical crack surfaces, the nodal releasing technique is used. In order to represent the material separation, the element connectivity at the current tip node must be modiﬁed; a duplicate node is created at the same location where the current tip resides. Few elements are released from the current tip node, and are assigned with a new duplicate node. This, in turn, generates new crack surfaces. In order not to destroy the FE mesh during this process, a consistent way of determining which elements getduplicate nodesisused. In this procedure, two angles are deﬁned: one is φ already deﬁned in Fig. 4, and the other is the angle formed by the negative normal at crack tip to the propagation vector (φ nn in Fig. 5a). Then, for each element, the angle (φ ) formed by the line connecting current tip to the centroid of the element and the normal is computed. The element is released and gets the duplicate node, if φ ∈ / [φ , φ ]. The elements that retain the current tip g p nn node are shaded in Fig. 5a. After nodal releasing and modifying element connectivity, the mesh close to the crack tip is plotted in Fig. 5b. At this point, the material separation is introduced and all the crack propagation operation are completed. Numerical examples Several examples of varying complexity are solved to demonstrate the eﬀectiveness of the proposed method. These examples exhibit single and mixed-mode behavior, involv- ing mono- and multimaterials. In order to closely examine the accuracy of the method, crack paths obtained from the present method are compared with experiments or results obtained from other methods in literature. The ﬁrst two examples consider stationary cracks, and the quantities calculated are compared with XFEM studies [30,31]. These examples consider highly nonlinear eﬀects evident from the crack tip blunting observed in the results. All the other examples involve complex crack propagation through the structure. Crack tip blunting Consider a single edge notched specimen with dimensions 2 mm × 6 mm. The crack occupies half-width as shown in Fig. 6a. The top surface is subjected to a ﬁxed displacement of 4 mm. All these details are taken from [30]. The strain energy function of the material is given by, = (tr C − 3) (15) NH The Lame parameters are set such that μ = 0.4225 MPa and the equivalent Poisson ratio in the linear regime ν = 0.49. When the material deforms, the crack surfaces move apart, and the initially sharp crack will blunt signiﬁcantly due to the material nonlinearity. The deformed conﬁguration of the structure is shown in Fig. 6b. The vertical displacement of crack surface nodes are plotted against their horizontal position in the reference state in Fig. 6c; for comparison, XFEM simulation results are taken from [30]. It is directly evident that the results obtained Sudhakar and Wall Adv. Model. and Simul. in Eng. Sci. (2017) 4:2 Page 13 of 23 coarse ﬁne 3.5 XFEM 2.5 1.5 0.5 -1 0 1 x(mm) a b c Fig. 6 Crack tip blunting. a Geometry. All dimensions are in mm. b Deformed conﬁguration. c Plot of vertical displacement of crack surface nodes against their horizontal position in reference state. XFEM results for comparison are taken from [30] from our simulations are matching well with the reported results. Moreover, the simu- lations using coarse and ﬁne mesh yield identical values, which shows that the reported results are converged with mesh density. The coarse and ﬁne mesh contain 1200 and 4800 uniform Cartesian elements respectively. It is to be mentioned that the XFEM study [30] was focused on incompressible materials, and the present results closely resembles the incompressible condition by taking ν = 0.49. J-integral computation J-integral is a crucial parameter in our work because both the crack propagation and crack kinking criterion are entirely based on this single parameter. In order to study the accuracy of J-integral evaluation, we consider an edge crack specimen under simple extension. A 2mm × 2 mm plate with μ = 0.4425 MPa and the equivalent Poisson ratio in the linear regime ν = 0.49 is taken, and the crack occupies half-width of the specimen. The strain energy function and boundary conditions are same as that of the previous example. These details are taken from an XFEM study [31], which reports the value of J-integral for large stretch ratios (λ). λ is deﬁned as the ratio of deformed length to the original length. We intentionally chose [31] as the reference for our validation because comparing J-integral values at high λ could be challenging. It can be seen from Fig. 7a that the J-integral values computed from our method are in excellent agreement with the XFEM results [31] even at large λ. The coarse and ﬁne mesh indicated in Fig. 7a contain 728 and 1600 uniform Cartesian elements respectively. At λ = 2.5, the diﬀerence between J-integral values computed using the coarse and ﬁne u (mm) y Sudhakar and Wall Adv. Model. and Simul. in Eng. Sci. (2017) 4:2 Page 14 of 23 coarse ﬁne Rashetnia, Mohammadi (2015) 1.5 0.5 1 1.5 2 2.5 1.3 coarse λ =2 ﬁne λ =2 coarse λ =2.5 1.2 ﬁne λ =2.5 1.1 0.9 0.8 0.7 23456 7 89 10 No. of layers Fig. 7 Computation of J-integral: a Variation with respect to λ. b Domain independency mesh is only 1.3%. At the same λ, the diﬀerence between the value reported in [31]and the present simulation using the ﬁne mesh is as low as 3.4%. This quantitative comparison shows that the J-integrals computed in our work are very accurate even at large λ. The variation of J with respect to the number of layers chosen around the crack tip as the integration domain is given in Fig. 7b. It can be seen that after initial changes, the value of J is stable after 5 layers, after which only minute variations exist in J.Inall the simulations presented in this work, 5 layers of elements around the crack tip are chosen to compute the J-integral. Having simulated a stationary crack, the remaining examples consider complex crack propagation through the structure. The comparison between the present results and the results obtained from the literature demonstrates the accuracy of the method. For all the following examples, the strain energy function is given by Eq. 4. Single edge cracked plate under mixed-mode loading In this example, as shown in Fig. 8a, the plate which is ﬁxed at the bottom edge, is subjected to shear stress τ = 1 on the top. The initial edge crack length is half of the plate width. The J-integral ref Sudhakar and Wall Adv. Model. and Simul. in Eng. Sci. (2017) 4:2 Page 15 of 23 ab Present Phongthanapanich, Dechaumphai (2004) Rao, Rahman (2000) 8.5 7.5 01 23 456 Fig. 8 Single edge crack plate under mixed mode loading: a Geometry, material parameters and loading conditions; all dimensions are in mm. b Contours of displacement in vertical direction. c Crack tip trajectory Lame parameters are set such that E = 30 MPa and ν = 0.25. The computational domain is discretized with 2736 elements with the whole area of crack propagation discretized with a ﬁne mesh. In this simulation, upon loading the crack propagates along a slightly curved path until it reaches the other end of the plate. In order to provide a detailed comparison, a zoomed y Sudhakar and Wall Adv. Model. and Simul. in Eng. Sci. (2017) 4:2 Page 16 of 23 view of crack path is plotted in Fig. 8c, and the result obtained from the present simulation is compared with two other studies: one utilizes a meshless method [53], and the other study is based on adaptive FEM [6]. It can be seen that the predicted crack tip trajectory matches very well with the results obtained from the other two studies. Crack in a drilled plate To demonstrate further the accuracy of the proposed method to simulate the crack path, the example given in [7] is considered. It reported the propagation of a crack from an initial notch in a beam which has three drilled holes. The study carried out both experimental and numerical tests, and observed a curvilinear crack propagation within the drilled plate. The geometrical conﬁguration, material properties, and the loading conditions are given in Fig. 9a.The Lame parameters are set such that E = 3GPa and ν = 0.35. In this example, the stress/strain ﬁelds are inﬂuenced by the presence of holes in the beam, and this provides interesting curvilinear crack tip trajectories. There are two simulation cases considered based on the location of the initial notch. These are dictated by the choice of a and b in Fig. 9a whose values are given in Table 1 for simulation-1 and simulation-2. As reported in [7], the crack path follows diﬀerent trajectories based on the choice of a and b, which are described as follows. Simulation-1 The location of the initial notch is given by a = 5and b = 1.5 mm. The crack is initially attracted towards the bottom hole, propagates near this hole, and got deﬂected away to end in the middle hole as shown in Fig. 9b. This is in accordance with the experimental results of [7], and other numerical studies [3,6,54]. Comparison with the experimental results show that the present simulation produces very good results; even the crack deﬂection near the bottom hole is predicted well in the simulation as can be directly seen from Fig. 9b. This is one of the very challenging validation test cases, owing to the complex crack tip trajectory involved. The developed methodology can be said to be accurate as it produces results that are matching very well with the experimental values even for this complex conﬁguration. Simulation-2 In this example, for which a = 6and b = 1 mm, the crack is attracted towards the middle hole, and directly ends in it (Fig. 9c). There are no crack deﬂections observed, and for this example as well, the results match excellently with the experiment (Fig. 9c). Four point beam with two notches In order to test the performance of the present method to simulate multiple cracks in a structure, the four point bending beam with two pre-existing notches, shown in Fig. 10a, is simulated. The beam is supported from below at two points, and is loaded at two other points. The material properties are also given in Fig. 10a. The computational domain is discretized with 7400 elements. This example is proposed by Bocca et. al. [55] who performed experiments on the structure, and also simulated them numerically. The crack paths through the FE mesh is given in Fig. 10b. To demonstrate the accuracy of the simulation, the crack paths obtained from the present method are compared with Sudhakar and Wall Adv. Model. and Simul. in Eng. Sci. (2017) 4:2 Page 17 of 23 Fig. 9 Bittencourt’s drilled plate problem. a Geometry, b Simulation-1, c Simulation-2. Experimental values for comparison are taken from [7] the results reported using a meshless method that incorporates crack tip singular ﬁelds as enrichments [56]; results presented for the ﬁnest meshless node distribution is used for the comparison. The comparison of crack paths is plotted in Fig. 10c. It can be seen that for both crack tips, the tip trajectory obtained from the present simulations matches very well with the reported value. Sudhakar and Wall Adv. Model. and Simul. in Eng. Sci. (2017) 4:2 Page 18 of 23 Table 1 Geometric parameters deﬁning notch location for Bittencourt’s drilled plate problem shown in Fig. 9a Simulation a b 15.0 1.5 26.0 1.0 0.2P P E=27GPa ν=0.18 80 P=-100N 320 80 400 Present Rabczuk, Zi (2006) 250 300 350 400 450 500 550 x(mm) Fig. 10 Four point bending beam with two notches: a Geometry, material parameters and loading conditions (not to scale). All dimensions are in mm, b crack path, c comparison of crack tip trajectories y(mm) 200 Sudhakar and Wall Adv. Model. and Simul. in Eng. Sci. (2017) 4:2 Page 19 of 23 Though our method can handle multiple cracks, care must be taken to make sure that the domains used to evaluate J-integral do not intersect with each other. However, this is not a speciﬁc drawback associated with our method alone. This is a common issue with other available approaches as well. Crack deﬂection due to inclusion Crack growth in the presence of an inclusion is studied in this example. Geometry, loading, and boundary conditions are given in Fig. 11a; they are taken from [57]. The conﬁguration consists of a rectangular plate which contains an oﬀ-centre circular inclusion. The Lame parameters of the plate are set such that E = 20 MPa and ν = 0.3. The objective of this plate study is to check whether the method is capable of accurately predicting the inﬂuence of this inclusion on crack propagation, which is already reported in [52,57]. The inclusion is characterized by the ratio of Young’s modulus of the plate to that of the inclusion (r = E /E ). Two values are considered; r = 10 which means that the plate incl. Young’s modulus of the inclusion is 10 times lower than that of the plate which is referred to as “soft” inclusion, and r = 0.1 that is referred to as “hard” inclusion. The Poisson ratio Fig. 11 Crack deﬂection due to inclusion a Geometry and loading conditions, b crack path Sudhakar and Wall Adv. Model. and Simul. in Eng. Sci. (2017) 4:2 Page 20 of 23 is assumed to be the same as that of the plate. The whole structural domain is discretized with 3213 elements. The eﬀect of the inclusion on the crack tip trajectory is shown in Fig. 11b. For soft inclusion, the crack is attracted towards the side of the inclusion; however, the crack does not end in it. In case of a hard inclusion, the crack deﬂects away from it. These observations are consistent with the already reported results [52,57]. Nonlinear elastic plate with a hole The above examples considered crack propagation with little material nonlinearity. This is evident from the fact that the crack remains sharp even after several propagation steps. The following example considers crack propagation involving high material nonlinearity under large deformation. A small oﬀ-centre hole is introduced in the geometric conﬁguration considered for the ﬁrst example, and this simulation allows the crack to propagate through the material. The Lame parameters are set to yield E = 10 GPa, ν = 0.3; critical J-integral, −2 J = 50 kJm . The top surface is subjected to a displacement of 0.5 mm. The geometric conﬁguration of this example is presented in Fig. 12a. As with the linear elastic examples, the loading (or the corresponding Dirichlet boundary condition here) is increased very smoothly from the zero initial value so that the inﬂuence of inertia is neglected. When the material starts deforming, as expected, the crack starts to blunt, and the J-integral value starts to increase. When J reaches J , then the crack starts to propagate; the deformed conﬁguration of the structure at which the crack starts propagating is depicted in Fig. 12b. Due to the presence of the hole, the crack slightly deﬂects upwards, as can be seen from Fig. 12b, c. From all these plots, one can infer that the crack tip is always blunt owing to the material nonlinearity, and the present method is able to model fracture behavior in such scenarios. Fig. 12 Plate with a hole. a Geometry (not to scale). All dimensions are in mm. b The conﬁguration at which the crack starts propagating. c An intermediate conﬁguration. d Final conﬁguration Sudhakar and Wall Adv. Model. and Simul. in Eng. Sci. (2017) 4:2 Page 21 of 23 Conclusion A ﬁnite element methodology to model mixed-mode crack propagation through non- linear elastic materials is proposed in this work. The striking feature of this method is that it facilitates, with minimal implementation eﬀorts, to update an existing large scale structural mechanics solver into a robust tool to handle single and multiple cracks. The method involves two steps: in the ﬁrst step, the governing equations of the structure are solved using nonlinear FEM by freezing the crack in the structure; in the next step, the solution obtained from the FEM is used to propagate the crack based on the maximum energy release rate criterion. Advancing the crack through a FE mesh requires a continual change in topology of the mesh, which is achieved in this work by utilizing a mesh reﬁt- ting approach. This method, as the name suggests, reﬁts the mesh at each instant of crack advancement in such a way that the crack propagates through an existing edge in the mod- iﬁed mesh. The mesh deformation strategies (for example used in ALE based methods) usually result in mesh tangling issues when attempting to handle topology changes in the mesh. This problem is circumvented in this work by splitting the quadrilateral elements into triangular elements in the crack tip neighborhood, which allows the possibility of local mesh connectivity to be modiﬁed. This step is crucial to preserve the quality of the mesh throughout the simulation, without which the mesh movement methods will fail. Examples involving single- and multi-materials with one or multiple cracks are reported. The obtained results are compared with experimental and other available computational methods. The comparison demonstrated that the present method accurately predicted the fracture behavior of all the examples considered. Authors’ contributions YS developed the idea, conducted numerical experiments and wrote draft. WAW ﬁne-tuned the research idea, suggested numerical experiments and revised the paper. Both authors read and approved the ﬁnal manuscript. Competing interests The authors declare that they have no competing interests. Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional aﬃliations. Received: 30 August 2016 Accepted: 5 June 2017 References 1. Sudhakar Y, Wall WA. A strongly coupled partitioned approach for ﬂuid-structure-fracture interaction. Int J Numer Methods Fluids; Submitted. 2. Wawrzynek PA, Ingraﬀea AR. An edge-based data structure for two-dimensional ﬁnite element analysis. Eng Comput. 1987;3:13–20. 3. Miehe C, Gürses E. A robust algorithm for conﬁgurational-force-driven brittle crack propagation with R-adaptive mesh alignment. Int J Numer Methods Eng. 2007;72:127–55. 4. Bouchard PO, Bay F, Chastel Y, Tovena I. Crack propagation modelling using an advanced remeshing technique. Comput Methods Appl Mech Eng. 2000;189:723–42. 5. Miranda ACO, Meggiolaro MA, Castro JTP, Martha LF, Bittencourt TN. Fatigue life and crack path predictions in generic 2D structural components. Eng Fract Mech. 2003;70:1259–79. 6. Phongthanapanich S, Dechaumphai P. Adaptive Delaunay triangulation with object-oriented programming for crack propagation analysis. Finite Elem Anal Design. 2004;40:1753–71. 7. Bittencourt TN, Wawrzynek PA, Ingraﬀea AR. Quasi-automatic simulation of crack propagation for 2D LEFM problems. Eng Fract Mech. 1996;55:321–34. 8. Camacho GT, Ortiz M. Computational modeling of impact damage in brittle materials. Int J Solids Struct. 1996;33:2899– 9. Ortiz M, Pandolﬁ A. Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation anal- ysis. Int J Numer Methods Eng. 1999;44:1267–82. 10. de Borst R. Numerical aspects of cohesive-zone models. Eng Fract Mech. 2003;70:1743–57. 11. Turon A, Dávila CG, Camanho PP, Costa J. An engineering solution for mesh size eﬀects in the simulation of delami- nation using cohesive zone models. Eng Fract Mech. 2007;74:1665–82. Sudhakar and Wall Adv. Model. and Simul. in Eng. Sci. (2017) 4:2 Page 22 of 23 12. Park K, Paulino GH. Cohesive Zone Models: a critical review of traction-separation relationships across fracture surfaces. Appl Mech Rev. 2013;64:060802. 13. Allix O, Ladevze P, Gilletta D, Ohayon R. A damage prediction method for composite structures. Int J Numer Methods Eng. 1989;27:271–83. 14. Genet M, Marcin L, Ladevze P. On structural computations until fracture based on an anisotropic and unilateral damage theory. Int J Damage Mech. 2014;23:483–506. 15. Celes W, Paulino GH, Espinha R. A compact adjacency-based topological data structure for ﬁnite element mesh representation. Int J Numer Methods Eng. 2005;64:1529–56. 16. Moës N, Dolbow J, Belytschko T. A ﬁnite element method for crack growth without remeshing. Int J Numer Methods Eng. 1999;46:131–50. 17. Sukumar N, Moës N, Moran B, Belytschko T. Extended ﬁnite element method for three-dimensional crack modelling. Int J Numer Methods Eng. 2000;48:1549–70. 18. Réthoré J, Gravouil A, Combescure A. An energy-conserving scheme for dynamic crack growth using the eXtended ﬁnite element method. Int J Numer Methods Eng. 2005;63:631–59. 19. Xiao QZ, Karihaloo BL. Improving the accuracy of XFEM crack tip ﬁelds using higher order quadrature and statistically admissible stress recovery. Int J Numer Methods Eng. 2006;66:1378–410. 20. Gupta V, Kim DJ, Duarte CA. Analysis and improvements of global-local enrichments for the generalized ﬁnite element method. Comput Methods Appl Mech Eng. 2012;245–256:47–62. 21. Gupta V, Duarte CA, Babuška I, Banerjee U. A stable and optimally convergent generalized FEM (SGFEM) for linear elastic fracture mechanics. Comput Methods Appl Mech Eng. 2013;266:23–39. 22. Nagarajan A, Mukherjee S. A mapping method for numerical evaluation of two-dimensional integrals with 1/r singularity. Comput Mech. 1993;12:19–26. 23. Béchet E, Minnebo H, Moës N, Burgardt B. Improved implementation and robustness study of the X-FEM for stress analysis around cracks. Int J Numer Methods Eng. 2005;64:1033–56. 24. Laborde P, Pommier J, Renard Y, Salaün M. High-order extended ﬁnite element method for cracked domains. Int J Numer Methods Eng. 2005;64:354–81. 25. Park K, Peraira JP, Duarte CA, Paulino GH. Integration of singular enrichment functions in the generalized/extended ﬁnite element method for three-dimensional problems. Int J Numer Methods Eng. 2009;78:1220–57. 26. Mousavi SE, Sukumar N. Generalized Duﬀy transformation for integrating vertex singularities. Comput Mech. 2010;45:127–40. 27. Mousavi SE, Sukumar N. Generalized Gaussian quadrature rules for discontinuities and crack singularities in the extended ﬁnite element method. Comput Methods Appl Mech Eng. 2010;199:3237–49. 28. Minnebo H. Three-dimensional integration strategies of singular functions introduced by the XFEM in the LEFM. Int J Numer Methods Eng. 2012;92:1117–38. 29. Dolbow J, Devan A. Enrichment of enhanced assumed strain approximations for representing strong discontinuities: addressing volumetric incompressibility and the discontinuous patch test. Int J Numer Methods Eng. 2004;59:47–67. 30. Legrain G, Moës N, Verron E. Stress analysis around crack tips in ﬁnite strain problems using the eXtended ﬁnite element method. Int J Numer Methods Eng. 2005;63:290–314. 31. Rashetnia R, Mohammadi S. Finite strain fracture analysis using the extended ﬁnite element method with new set of enrichment functions. Int J Numer Methods Eng. 2015;102:1316–51. 32. Donea J, Huerta A, Ponthot JP, Rodrìguez-Ferran A. Arbitrary Lagrangian-Eulerian methods. In: Stein E, Borst R, Hughes TJR, editors. Encyclopedia of computational mechanics, vol. 1. Hoboken: Wiley; 2004. 33. Koh HM, Haber RB. A mixed Eulerian-Lagrangian model for the analysis of dynamic fracture. Illinois: Dept. of Civil Engineering, University of Illinois at Urbana-Champaign, Civil Engineering Studies, Structural Research Series; 1986. p. 34. Koh HM, Lee HS, Haber RB. Dynamic crack propagation analysis using Eulerian-Lagrangian kinematic descriptions. Comput Mech. 1988;3:141–55. 35. Abdelgalil AI. Modeling of dynamic fracture problems using ALE ﬁnite element formulation. Vancouver: The University of British Columbia; 2002. http://hdl.handle.net/2429/13320. 36. Amini MR, Shahani AR. Finite element simulation of dynamic crack propagation process using an arbitrary Lagrangian Eulerian formulation. Fatigue Fract Eng Mater Struct. 2013;36:533–47. 37. Bruno D, Greco F, Lonetti P. A fracture-ALE formulation to predict dynamic debonding in FRP strengthened concrete beams. Compos Part B. 2013;46:46–60. 38. Ponthot JP, Belytschko T. Arbitrary Lagrangian-Eulerian formulation for element-free Galerkin method. Comput Meth- ods Appl Mech Eng. 1998;152:19–46. 39. Rashid MM. The arbitrary local mesh replacement method: an alternative to remeshing for crack propagation analysis. Comput Methods Appl Mech Eng. 1998;154:133–50. 40. Tabiei A, Wu J. Development of the DYNA3D simulation code with automated fracture procedure for brick elements. Int J Numer Methods Eng. 2003;57:1979–2006. 41. Browne PA, Budd CJ, Piccolo C, Cullen M. Fast three dimensional r-adaptive mesh redistribution. J Comput Phys. 2014;275:174–96. 42. Budd CJ, Russell RD, Walsh E. The geometry of r-adaptive meshes generated using optimal transport methods. J Comput Phys. 2015;282:113–37. 43. Belytschko T, Liu WK, Moran B. Nonlinear ﬁnite elements for continua and structures. Hoboken: Wiley; 2001. 44. Bonet J, Wood RD. Nonlinear continuum mechanics for ﬁnite element analysis. Cambridge: Cambridge University Press; 1997. 45. Rice JR. A path independent integral and the approximate analysis of strain concentration by notches and cracks. J Appl Mech. 1968;35:379–86. 46. Moran B, Shih CF. Crack tip and associated domain integrals from momentum and energy balance. Eng Fract Mech. 1987;27:615–42. Sudhakar and Wall Adv. Model. and Simul. in Eng. Sci. (2017) 4:2 Page 23 of 23 47. Shivakumar KN, Raju IS. An equivalent domain integral method for three-dimensional mixed-mode fracture problems. Eng Fract Mech. 1992;42:935–59. 48. Ma L, Korsunsky AM. On the use of vector J-integral in crack growth criteria for brittle solids. Int J Fract. 2005;133:L39–46. 49. Erdogan F, Sih GC. On the crack extension in plates under plane loading and transverse shear. J Basic Eng. 1963;85:519– 50. Sih GC. Strain-energy-density factor applied to mixed-mode fracture problems. Int J Fract. 1974;10:305–21. 51. Hussain MA, Pu SL, Underwood JH. Strain energy release rate for a crack under combined mode I and mode II. Fract Anal ASTM STP. 1974;560:2–28. 52. Bouchard PO, Bay F, Chastel Y. Numerical modelling of crack propagation: automatic remeshing and comparison of diﬀerent criteria. Comput Methods Appl Mech Eng. 2003;192:3887–908. 53. Rao BN, Rahman S. An eﬃcient meshless method for fracture analysis of cracks. Comput Mech. 2000;26:398–408. 54. Areias P, Dias-da-Costa D, Alfaiate J, Júlio E. Arbitrary bi-dimensional ﬁnite strain cohesive crack propagation. Comput Mech. 2009;45:61–75. 55. Bocca P, Carpinteri A, Valente S. Mixed mode fracture of concrete. Int J Solids Struct. 1991;27:1139–53. 56. Rabczuk T, Zi G. A meshfree method based on the local partition of unity for cohesive cracks. Comput Mech. 2006;39:743–60. 57. Natarajan S. Enriched ﬁnite element methods: advances & applications. Cardiﬀ: Cardiﬀ University; 2011.
"Advanced Modeling and Simulation in Engineering Sciences" – Springer Journals
Published: Dec 1, 2017
You can share this free article with as many people as you like with the url below! We hope you enjoy this feature!
Read and print from thousands of top scholarly journals.
Already have an account? Log in
Bookmark this article. You can see your Bookmarks on your DeepDyve Library.
To save an article, log in first, or sign up for a DeepDyve account if you don’t already have one.
Copy and paste the desired citation format or use the link below to download a file formatted for EndNote
Access the full text.
Sign up today, get DeepDyve free for 14 days.
All DeepDyve websites use cookies to improve your online experience. They were placed on your computer when you launched this website. You can change your cookie settings through your browser.