Dynamic crack growth based on moving mesh method

Composites Part B 174 (2019) 107053

Contents lists available at ScienceDirect

Composites Part B journal homepage: www.elsevier.com/locate/compositesb

Dynamic crack growth based on moving mesh method Francesco Fabbrocino a, Marco Francesco Funari b, Fabrizio Greco b, *, Paolo Lonetti b, Raimondo Luciano c, Rosa Penna d a

Department of Civil Engineering, Pegaso University, Naples, Italy Department of Civil Engineering, University of Calabria, Rende, Italy c Department of Engineering, University of Naples Parthenope, Naples, Italy d Department of Civil Engineering, University of Salerno, Salerno, Italy b

A R T I C L E I N F O

A B S T R A C T

Keywords: Crack propagation ALE Moving mesh technique Finite element method Composite structures

A new methodology to predict dynamic crack propagation under generalized loading conditions is proposed. The numerical modeling combines structural mechanics and moving mesh method with the purpose to predict ge­ ometry variation produced by the evolution of existing material discontinuities. In particular, moving mesh method is implemented to enforce crack tip displacements by using an explicit crack criterion based on refer­ ential and moving configurations. In this framework, the use of mesh regularization method based on proper rezoning equations is able to reduce the use of remeshing attempts, typically required by standard crack prop­ agation procedures. Dynamic crack growth is predicted by a rate dependent criterion, expressed in terms of crack angle and driven forces based on energy release rate definition. The model is quite suitable to predict the evolution of material discontinuities, typically observed in composite structures. Numerical implementation, developed in the framework of a finite element formulation and details on the solving procedure, are presented. The proposed modeling is validated by several comparisons with experimental and numerical data, which show accuracy and robustness of the numerical approach. Moreover, sensitivity analyses in terms of mesh dependence and time required for the solving procedure are also developed.

1. Introduction Dynamic fracture mechanics is considered an active area of research, since the simulation crack growth phenomena affects many fields, ranging from structural or mechanical engineering, earthquake wave propagation and high speed impact/contact phenomena. Typically, the failure of brittle materials, such as composite structures, is dominated by the formation or evolution of internal material discontinuities, whose simulation requires advanced techniques to predict high crack tip speeds or complex phenomena such as branching or coalescence mechanisms [1]. From macroscopic point of view, the structural system is affected by a considerable loss of stiffness, which may lead the structure to the complete failure, catastrophically. In order to simulate crack propaga­ tion, the field of Computational Fracture Mechanics (CFM), specifically in dynamics, is currently object of several studies and formulations. Most of the models available from the literature are developed in the framework of the Finite Element Method (FEM), because of its versatility and simplicity to model complex structures, approximating the solutions

of partial differential equations [2]. Classical FE models are based on the definition of local constitutive approaches, mainly devoted to identify the stiffness degradation produced by the crack growth, such as for instance Continuum Damage Mechanics (CDM) [3]. However, since such models are affected by numerical instabilities produced, essen­ tially, by their ill-posed boundary value problems, numerical strategies based on high order continuum models and non local approaches are proposed [4]. Contrarily to previous models, FEM is also utilized to simulate the behavior of explicit material discontinuities in terms of crack propagation and resistance curves. To this end several numerical methodologies are proposed, in which the crack is predicted by releasing existing nodes on specific boundaries, where the crack is expected to advance [5] or introducing proper stiffness penalization factors in the constitutive laws to eliminate FEs affected by the crack path [6]. Such approaches require a very fine mesh around the moving crack tip with a notable increase of the computational costs of the model. Alternatively, Cohesive Zone Models (CZM) are frequently utilized in the literature with large success especially in those applications in which the crack is

* Corresponding author. Department of Civil Engineering University of Calabria, 87036, Rende, Cosenza, Italy. E-mail addresses: [email protected] (F. Greco), [email protected] (P. Lonetti). https://doi.org/10.1016/j.compositesb.2019.107053 Received 7 May 2019; Received in revised form 17 June 2019; Accepted 18 June 2019 Available online 22 June 2019 1359-8368/© 2019 Elsevier Ltd. All rights reserved.

F. Fabbrocino et al.

Composites Part B 174 (2019) 107053

known a priori, i.e. where the path coincides with an external boundary [7]. A generalization of such procedures is achieved by introducing the interfaces at interelement boundaries, in which an arbitrary crack path is able to grow in a self consistent way. However, CZMs are affected by mesh dependence problems and numerical instabilities due to the soft­ ening constitutive laws, produced by the presence of a high number of elements involved in the mesh grid points. Moreover, in dynamics, the presence of an initial stiffness of the cohesive elements affects internal wave speeds leading to spurious oscillation phenomena of the traction forces [8]. Problems reported to above can be circumvented by means of adaptive mesh procedures, which modify by mesh connectivity to reproduce the crack evolution [9]. The major problems of such ap­ proaches are related to mesh connectivity constraints and topology consistency to be achieved during the crack growth [10]. The mesh should be modified by generation algorithms, which verify the interplay among mesh elements in the structure. Additional complexities are observed in those cases in which multilevel adaptive meshes are implemented, which decrease computational accuracy, numerical convergence and efficiency [11,12]. Previous models are developed in the framework of microscale, in which an explicit crack geometric representation is simulated. In order to reduce mesh dependence problems, related to the accuracy of the crack tip discretization, kinematic nongeometrical approaches, such as Extended or Generalized Finite Element Methods (XFEM, GFEM) are proposed. In this framework, the crack is simulated by enriching strain or displacement functions to include internal discontinuities produced by the presence of microcracks [13,14]. Therefore, the initial defects are incorporated in the model at arbitrary points, in which, potentially, the crack can be triggered without defining in advance a crack path. How­ ever, in order to ensure accuracy and efficiency in the prediction of fracture variables, a refined mesh is required in all the elements affected, potentially, by the presence of microcracks. This may lead global mesh enrichment, which inevitably produces expensive computational modeling. In addition, the crack path is based on specific geometry mesh intersection criteria, which require the implementation additional pro­ gramming and calculation procedures to standard FE software. A hybrid approach with respect to geometrical or nongeometrical representations is defined by using moving mesh methodologies based on Arbitrary Lagrangian-Eulerian (ALE) formulation. Typically, appli­ cations are implemented in the framework of fluid mechanics to reduce mesh distortions [15,16] or to take into account moving boundaries in Fluid-Structure interaction (FSI). However, a generalization in the framework of solid mechanism and subsequently in Fracture Mechanics was also developed in the literature. In particular, the basic idea is to reproduce crack motion as geometric variation of the interface regions, by means of specific rezoning or regularization methods [17–19]. Despite existing numerical methods utilized in CFM, such methodology appears to handle more efficiently the crack growth, since it reduces strongly the application of remeshing procedures. However, specific additional tasks should be accomplished to apply this numerical strategy to complex structures. In particular, most of the existing models in the literature are developed under in static loading conditions, neglecting a priori any initial effect triggered by the moving crack tip. As shown in Refs. [20,21], such assumption may produce notable underestimations of the crack growth variables with respect to experimental observations. Moreover, the crack path is assumed to be one dimensional in the case of interfacial crack growth or, in 2D/3D cases with sharp shape. The main goal of this paper is to generalize the numerical implementation pro­ posed in Ref. [22] by using proper time dependent crack criteria. This allows the development of a series of numerical results demonstrating the effectiveness of the method to simulate crack growth in continuum media under dynamic loading conditions or impact phenomena. The outline of the paper is as follows. Section 2 presents an overview of ALE formulation, whereas in Section 3 the governing equations for the ALE and Fracture Mechanics problems are presented. Finally, in Section 4, comparisons with dynamic and impact loading schemes are reported to

validate the proposed modeling. 2. Overview of the ALE formulation A two dimensional domain Ω and the corresponding external contour

∂Ω, in which prescribed displacements and traction conditions on ∂Ωu and ∂Ωt , respectively, are considered. Moreover, Γ is an internal mate­

rial discontinuity, which potentially may grow along the normal cone outgoing from the tip front direction. The model is based on a linear elastic isotropic material in small deformations. Consistently with ALE approach, two different configurations are supposed, defined as Refer­ ential (R) and Moving (M) ones (Fig. 1), which identify, for each mesh point, the mapping between the current (moving) and fixed (referential) nodes: (1)

X ¼ ΦðX ; tÞ X ¼ Φ 1 ðX ; tÞ eM e eR eR e eM

� � where Φ : CR →CM with ΦT ¼ ΦX ΦY is assumed to be invertible e e with continuous inverse. The spatial derivative of a generic function f ¼ e f ðX ; tÞ is defined on the basis of the Jacobian matrix of the trans­ e eM formation (J) and it is expressed by means of the following relationships: e 3 2 ∂ΦX ∂ΦY 7 6 6 ∂XR ∂XR 7 7 r f ¼ J 1 ⋅r f with J ¼ 6 (2) 6 7 e 4 ∂ΦX ∂ΦY 5 e Xe e e Re ∂YR ∂YR In addition material and referential derivatives are defined as the rate of change of the function with X and X fixed, respectively: eM eR d f_ ¼ f ðX ; tÞjX ; dt e M e

0

f ¼

M

d f ðX ; tÞj dt e R eX

(3)

R

where the dot and the prime are used here and in the following to represent material and referential time rates, respectively. According to ALE formulation, the time derivatives of a generic physical field, in referential and material configurations can be related by the following relationships: 0 f_ ¼ f

X

0

M

d d 0 f ðX ; tÞ with X M ¼ ΦðX ; tÞjX dX e e M dt e e R e R eM

(4)

where X M represents the relative velocity of the grid points in the ma­ terial reference system. Analogously, starting from Eq. (4), the second time derivative is evaluated recursively as: 0

f€¼ f 00

0

2r f ⋅X eX

0

M

r f X þ r ðr f ÞX e Xee M e X e Xe

0

M

X

0

M

þ r fr X e Xe e X

0

M

X

0

M

(5)

3. Governing equations The governing equations are derived on the basis of the generalized principle of virtual works: δU þ δK ¼ δV

(6)

where δU, δK and δV are virtual works of internal, inertial and external forces, respectively. In particular, previous quantities are defined on the current material configuration, but are written by the use of Eqs. (1)–(5) in the referential domain. In particular, the terms reported in Eq. (6) can be defined as follows:

2

F. Fabbrocino et al.

Composites Part B 174 (2019) 107053

Fig. 1. Referential (CR) and Moving (CM) configurations: relationship between X and X . eR eM

Z

Z δU ¼ h Ω

Eε⋅δεdA ¼ h ee e

ΩR

� � E J 1 r u ⋅ J 1 r u JdA e e e Re e e Re

(ERR) and propagation angle θ. In order to evaluate the crack growth phenomenon, an incremental-iterative procedure should be imple­ mented, in which the unknown quantities at the current time step are represented by the crack tip displacements and crack propagation angle. To this end, a fracture function depending from fracture tip variables, such as for instance ERR or SIF, is required. In addition, angle prediction should be evaluated since the crack growth is supposed to advance with arbitrary direction in a 2D region. Therefore, the crack growth is ob­ tained by solving the following equations:

(7)

Z �d2 u� e ⋅ δudA ¼ δK ¼ h ρ dt2 e Ω

Z

ρ u

¼h

e

ΩR

2J 1 r u⋅X e Re e M

� r uX þ J 1 r J 1 r u X X e R ee M eR e Re eM eM

(8)

� þ J 1 r J 1 r X X ⋅δuJdA eR e R eM eM e Z δV ¼ h ∂Ω

∂fG with λ_ F � 0 and λ_ F fG ¼ λ_ G fG ¼ 0 fG ðG; θ; G0 Þ � 0 θ ¼ gF ðGÞ δ_ F ¼ λ_ F ∂G (11)

Z q⋅δudS ¼ h q⋅δuJ s dS e e e e ∂ΩR

(9)

where λ_ F is the fracture multiplier, G represents the fracture variable, i.e. total ERR, θ is the angle of propagation depending from the crack kinking function gF and δ_ F is the crack tip displacement. It is worth noting that, fG or gF are not specified at this stage, since, in the literature, many crack growth or propagation angle criteria are developed depending from material characteristics and phenomenological aspects. The proposed model is not strictly dependent, specifically, from the crack growth or kinking criterion, since it is quite general to include several analytical expressions depending from the material and crack growth characteristics. Specific choices will be described in the valida­ tion procedure presented in the following section for comparison purposes. In order to describe the shape variation produced by the cracked surface, specific boundary conditions and field equations are required.

where h is the thickness of the structure, E is the classical elastic matrix, e � � u with uT ¼ U1 e e " #

εx εy

γ xy

U2

is the displacement vector, ε with εT ¼ e e

is the strain vector, ρ is the per unit area material den­

sity, ðJ; Js Þ are determinants of the area and boundary transformation matrixes from material to referential configurations and q is the traction e force vector applied on the loaded boundaries. The resulting equation, concerning the structural formulation, is defined by substituting Eq. (7)(9) into Eq. (6), leading to the following relationship:

Z

� � 2J 1 r u⋅X r uX þ J 1 r J 1 r u X X þ J 1 r J 1 r X X ⋅δuJdAþ e e R e e M e R ee M eR e Re e M e M eR e R eM eM e Ω ZR Z � � 1 1 E J r u ⋅ J r δu JdA ¼ q⋅ δuJ S dS e e e Re e e R e e e

ρ u

ΩR

(10)

∂ΩR

Introducing the nodal mesh displacement vector, defined as ΔX ¼ X e eM X , boundary conditions on the crack tip contour enforce normal crack eR tip extension according to Eq. (11), as follows: � � ΔX ⋅ n ¼ δ_ F with nT ¼ sinðθÞ cosðθÞ on SF (12) e e e

It is worth noting that Eq. (10) contains not only variables related to structural mechanics, but also the positions and corresponding time and spatial derivatives of mesh points, i.e X and X ; X00 M . Therefore, eM eM additional equations are required to enforce mesh motion on the basis of a crack growth criterion defined in terms of fracture variables. In particular, geometry variation is based on explicit kinematic boundary conditions on the external contours, where the mesh motion is fixed (Fig. 2). Moreover, the crack tip motion is defined on the basis of a growth criterion expressed as a function of the Energy Release Rate

where n is the direction vector at generic point of the cracked surface, SF e 3

F. Fabbrocino et al.

Composites Part B 174 (2019) 107053

Fig. 2. Schematic representation of the crack tip motion between the i-th and the i-thþ1 time step.

is a small region enclosing the crack tip. In addition, rezoning equations as well as proper boundary conditions are required to redistribute the material points adjoining the cracked surface. Moreover, in order to ensure accuracy in the prediction of fracture variables due to mesh distortions during the crack growth, the region SF is moved rigidly (Fig. 2). This is achieved by means of LMM, which is utilized to enforce the motion of the crack region by means of implicit boundary condi­ tions. The mesh evolution in the remaining part of the structure, not affected by internal boundary conditions, is described by means of regularization equations according to Laplace smoothing method. It is worth noting that the choice of this method is purely arbitrary from the theoretical point of view, since the formulation is quite general to comply other well-known regularization techniques. Therefore, the variational form of the governing equations of the ALE problem are defined by means of the following expressions: Z Z Z rX : rδX dV ðΔX n δ_ F ÞδΛdA rδΔX ΛdA ¼ 0 (13) e ee e ee e e e ee VR

SF

discretization matrices. The discretized equations are solved by using an implicit time integration scheme based on a variable time step-size and backward differentiation formula (BDF), implemented by using an external subroutine that interacts with a COMSOL MULTIPHYSICS FE software [23]. In Fig. 3, the computational procedure is presented, synoptically. In particular, at a generic time step, the unknown variables are both the displacements vector and the position of the mesh points. It is worth noting that the coupled procedure with respect to the operator splitting technique reproduces more accurately the advancing phe­ nomena. Moreover, during the time integration, owing to the fast speeds involved in the crack advance, a small time step size is utilized, which typically, no more than 1/10 of the fundamental period of vibration of the structure. In order to compute accurately the fracture variables, a fine discretization mesh and a standard numerical integration method (quadrature) is adopted on the contour line and on the area surrounding the crack tip. Since mesh distortions are introduced due to crack tip motion, the numerical model makes the use of adaptive remeshing formulation as far as lower mesh quality parameter reaches the corre­ sponding tolerance value. However, the recourse to adaptive remeshing procedure is quite limited in view of the ALE approach, which utilizes mesh regularization method to reduce mesh distortions during the crack tip motion. The proposed model is quite general and can be adopted to simulate quasi-static crack evolutions. A detailed description of the proposed formulation to solve static problem is reported in Ref. [22]. Such task is achieved by solving classical governing equations within a time inde­ pendent crack growth criterion.

SF

where Λ is the LMM vector. Eq. (10) and Eq. (13) with conditions defined e by Eq. (11) are implemented in a FE formulation, by using a proper discretization of the geometric domain into Ne elements and the use of isoparametric elements defined by the following set of structural and ALE variables: Ne X

Ne Ne Ne X X X N U; u¼ N U; u¼ N U; λ¼ N Λ I¼1 e I e I e I¼1 e I e I e I¼1 e I e I e I¼1 e λI e I N N N e e e X X X X ¼ N X ; X ¼ N X ;X ¼ N X ; eM e e e e e e MI M MI M MI MI I¼1 I¼1 I¼1 e MI e MI

u¼ e

(14)

4. Validation and numerical results In order to validate the proposed formulation, comparisons with numerical and experimental data under dynamic loading conditions are reported. In particular, results are proposed with the purpose to verify the prediction of proposed model in terms of mesh dependence, computational efficiency, and numerical complexity. At first, the analysis investigates the behavior of a polymethyl methacrylate (PMMA) rectangular plate subject to in-plane opening prescribed displacement at the upper and lower edges, which was analyzed both numerically and experimentally in Ref. [24]. The loading, boundary conditions and geometry configuration are represented in Fig. 4, whereas the mechanical parameters are summarized in Table 1. Moreover, the comparison is developed under a “application phase simulation” according to definition provided in Ref. [25], in which the crack propagation is simulated by using a relationship between fracture toughness and crack tip speed determined by means experimental data. In the present analysis, the crack function fG is expressed as follows:

where Ne represents the number of elements utilized, ðN; N ; N Þ are the e eM eλ matrix interpolation functions and ðU ; U ; U ; X ; X ; X ΛI Þ are the e I e I e I eM eM eM nodal vector functions. Substituting Eq. (14) into Eq.s (10) and Eq. (13), the following discrete equations defined by the summation over the ne elements, are obtained: Ne X k¼1 Ne X k¼1

M i1 U i þ e k ek W i ΔX i þ e k ek

Ne X k¼1

Ne X k¼1

0i

Ci U k þ ek

Ne X k¼1

Ki Ui þ ek e k

Qmi Λmi ¼ 0; ðALEÞ e k ek

Ne X k¼1

Pi ¼ 0; ðSTRUÞ ek

(15)

(16)

ðM Þ and ðK Þ are the consistent mass and stiffness matrices of the ek ek structural system including nonstandard terms due to positional vari­ ables, P is the external load vector, W and Qm are the ALE and LLM ei e e 4

F. Fabbrocino et al.

Composites Part B 174 (2019) 107053

Fig. 3. Schematic representation of the algorithm for dynamic crack propagation.

Fig. 4. PMMA rectangular plate: geometry, loading and intial mesh configuation adopted for the simulations. Table 1 PMMA rectangular plate: mechanical parameters. E ½MPa�

ν

μ ½MPa�

ρ ½kg m 3 �

G0c ½N m

3090

0.35

1144

1180

300

fG ¼

G Gc

where ðG; Gc Þ correspond to the computer ERR and related fracture toughness, respectively. According to Ref. [24] fracture toughness is expressed in terms of the crack tip speed by means of different re­ lationships for Rate Independent (RI) or Rate Dependent (RD) laws, which are defined by the following expressions: Gc ðcÞ ¼ � 1

G0c

� Rate Independent ðRIÞ

(18)

c cR

� Gc ðcÞ ¼ G0 log

cL

cL

� c

Rate Dependent ðRDÞ

(19)

where ðG0c ; cR Þ are material parameters corresponding to static fracture toughness and Rayleigh wave speed and ðG0 ; cL Þ are material initiation toughness and limit crack tip speed. In order to compute numerically the ERR, path independent J integral is developed consistently to the expression provided in Ref. [26], here reported for completeness: Z Z Jk ¼ h½ðW þ KÞnk ti ui;k �dS þ h ½ðρu€ fi Þui;k ρu_i u_i;k �dA (20) ∂SF

�

cr ½m s 906

1

�

G0 ½N m 2000

1

�

cL ½m s

1

�

675

the fracture variables. It is worth noting that material parameters introduced in Eq.18 and 19, reported in detail in Table 1, are taken from the literature to verify the capability of the proposed model, but their intrinsic definition is out of scope of the present analysis. Numerical and experimental data are collected by means of a two-step analysis. At first, the crack is taken as fixed and the structure is subjected to prescribed opening displacements (Fig. 1). At this stage, the total elastic energy of the structure is accumulated; once it reaches the prescribed value, the crack tip is free to grow. From the numerical point of view, a restart procedure from the initial elastic configuration is achieved and the crack growth is enforced by the activation of the crack growth criterion. Experimentally, the crack growth control is defined by the use of a razor cutter, which is applied to the specimen to produce crack advance. Several loading cases involving in different ranges of prescribed dis­ placements at the upper and lower edges are analyzed (Fig. 4). It is worth noting that, in the present loading scheme, the crack advances in a pure mode I, because the geometric and loading configuration of the structure. As a consequence, the crack path coincides with the sym­ metric horizontal axis, leading the crack angle always equal to zero. For each value of the prescribed displacement, the stored energy W0 is computed by using the following equation:

(17)

1�0

1

SF

W0 ¼

whereui ; u_i and u€i are the displacement velocity and acceleration of the material point, ρ is the material density, SF and ∂SF are surface and related contour enclosing crack tip, respectively. As shown in Fig. 2, SF represents the moving region (Fig. 2), which is able to follow rigidly the crack propagation ensuring a high level of accuracy in the prediction of

2EΔ20 H

(21)

where Δ0 is the prescribed edges displacement and H is the height of the rectangular plate. The investigation is performed with respect to the following loading configurations:

5

F. Fabbrocino et al.

� � � � � �

Δ0 ¼ 0:06 mm Δ0 ¼ 0:08 mm Δ0 ¼ 0:10 mm Δ0 ¼ 0:12 mm Δ0 ¼ 0:14 mm Δ0 ¼ 0:16 mm

Composites Part B 174 (2019) 107053

→ → → → → →

W0 W0 W0 W0 W0 W0

¼ 1390 N m ¼ 2470 N m ¼ 3860 N m ¼ 5562 N m ¼ 7570 N m ¼ 9890 N m

1 1 1 1 1 1

(L1) (L2) (L3) (L4) (L5) (L6)

The initial mesh configuration used for the analyses is generally coarse in the region externally located to the crack tip region. As shown in Fig. 4b, the tip region is meshed by using triangular plane stress el­ ements with characteristic length equal to ΔD =R ¼ 1=4 featuring Lagrange quadratic interpolation functions, whereas the remaining domain presents a mesh transition with maximum discretization length equal to ΔD =R ¼ 1=12, involving in total 5919 DOFs. The mesh gen­ eration procedure is consistent with a Delaunay tessellation method, which is able to minimize the interpolation error among all other tri­ angulations of the same set of vertices. In order to achieve accuracy in the time integration a maximum time step equal to 5E-7 is adopted. At first, RI analyses are discussed in Fig. 5, in which comparisons with numerical results are reported in terms of the crack tip speeds versus crack tip position. In the same figures, synoptic representations of the crack growth configurations are reported. From these analyses, as expected by the crack growth function, the crack tip speed increases with the stored energy, since a large amount of energy is available to trigger crack tip motion. The curves present largest values during the initiation phase, because of the instability nature of the dynamic pro­ cess. Subsequently, as far as the crack tip position increases, a steady state configuration with constant speed is reached. In Fig. 6, comparisons in terms of accumulated strain energy and crack tip speed are reported. All the points, including the one related to numerical data obtained in Ref. [24], show significant discrepancies with respect to the experimental results. In order to verify the capability of the proposed model and to better emphasize the influence of dynamic contributions included into the definition of the ERR, i.e. Eq. (20), two different solutions are reported, i.e. RI-st and RI-dy. In particular, the difference relies in the definition of the ERR, where only the static part is considered (RI-st) or the complete dynamic definition (RI-dy). The comparisons show that RI-st solution matches good with the numerical data obtained in Ref. [24]. However, RI-dy shows a better prediction, because of the kinetic contributions. However, the discrepancies be­ tween numerical and experimental data remain, but they can be

Fig. 6. RI law: Strain energy level vs average crack tip speed (W0-c/cr), com­ parisons with numerical and experimental data arising from Ref. [24].

Fig. 7. RD law: crack tip speed as a function of the crack tip position (c/cr-XT/ L), comparison with numerical data arising from Ref. [24].

explained by the presence of local micro-cracking generation produced at large crack speeds due to micro branching phenomena, which are responsible of an unpredicted fracture toughness of the material and energy dissipation. Previous results are extended and generalized in the case of a RD fracture criterion law. At first, the analysis is reported in terms of normalized values of crack tip speeds as a function of crack tip positions. Comparisons with numerical data, available from Ref. [24], are reported in Fig. 7, for several values of prescribed initial opening displacements. The results show a similar trend to the one reported in Fig. 5, with large speeds during the initiation phase and a regime steady state phase at large crack tip displacements. Globally, the speeds are quite smaller than those obtained in Fig. 5 and are always lower than the parameter speed cL. Fig. 8 analogously to previous case shows the comparisons in terms of W0 c points by introducing the average crack tip speeds obtained from the RD simulations. In particular, the results denote how, by adopting

Fig. 5. RI law: crack tip speed as a function of the crack tip position (c/cr-XT/L), comparison with numerical data arising from Ref. [24]. 6

F. Fabbrocino et al.

Composites Part B 174 (2019) 107053

Fig. 8. RD law: Strain energy level vs average crack tip speed (W0-c/cr), comparisons with numerical and experimental data arising from Ref. [24]. Fig. 10. RD law, loading configuation L3: comparisons in term CPU time and number of degrees of freedom for the numerical simulations performed by using different mesh discretizations (M1, M2 and M3).

the RD fracture toughness law, the proposed formulation matches perfectly both experimental and the numerical data reported in Ref. [24]. Such comparisons prove the accuracy of the RD approach, although it is strictly dependent by internal parameters such as ðG0 ; cL Þ, introduced in the crack function, which requires specific experimental procedures to be detected. Contrarily, the crack growth function utilized in RI comparisons, is only dependent by the static toughness and Ray­ leigh wave speed, which can be considered as known quantities depending from the material characteristics. Finally, in order to verify the computational efficiency of the pro­ posed model, a parametric analysis in terms of mesh discretization (L3) is developed. To this end the following mesh discretizations defined in terms mesh size ðΔDÞ and radius of Sf region (R), are considered: � characteristic length equal to ΔD =R ¼ 1=2 in the transition mesh in the remaining part of the plate length equal to ΔD =R ¼ 13=1 (M1). � characteristic length equal to ΔD =R ¼ 1=4 in the transition mesh in the remaining part of the plate length equal to ΔD =R ¼ 12=1 (M2). � characteristic length equal to ΔD =R ¼ 1=4 in the transition mesh in the remaining part of the plate length equal to ΔD =R ¼ 3=1(M3).

adopting a relatively coarse discretization, i.e. M1 and M2, it allows to save a notable amount of the computational time equal to 64% and 42% with respect to the M3 refined configuration. In Fig. 11, comparisons in terms of crack tip speed as a function of the crack tip position, for each mesh configurations (M1, M2, and M3), are reported. The proposed model is able to describe efficiently dynamic crack propagation phenomena by using a relatively coarse discretiza­ tion, i.e. M1. In any case, no instability or divergence phenomena are observed, since the fracture solution is stable and coincident with the refined one. Therefore, the numerical results are not affected by mesh dependent phenomena, since the mesh discretization in the moving domain (SF) region is able always to detect accurately the fracture var­ iables. Fig. 12 shows a synoptic representation of the mesh motion during the crack propagation for the M1 configuration. The result shows how the mesh movements of the crack tip region are enforced rigidly ensuring both the required accuracy in the prediction of the fracture variables and a reduced level of mesh dependency. According to ALE model, the elements of the Ω region are stretched due to the rezoning or regularization rules, generating a consistent discretization of the tran­ sition mesh and a strong reduction of the computational cost of the analysis. Additional comparisons are proposed for a complex case regarding an impact simulation in a PMMA specimen with a hole. This case study has been experimentally investigated by Ref. [27], and results were numerically adopted to verify the computational consistency of several

tip region and with maximum tip region and with maximum tip region and with maximum

Fig. 9 shows a zoomed view of the mesh discretization around the crack tip region for each numerical grid points adopted in the analysis (M1, M2, and M3). Fig. 10 investigates the computational efficiency in terms of both CPU time and number of DOFs. The results show how

Fig. 9. Mesh discretization detail around the crack tip for M1, M2 and M3 configurations. 7

F. Fabbrocino et al.

Composites Part B 174 (2019) 107053

assumed to be known, whereas crack direction is evaluated at current crack tip displacement by means of a crack angle criterion. In the present analysis, the crack speed is explicitly prescribed by adopting a piecewise function, reported in Fig. 15, while the crack is assumed to propagate along the direction computed by using the crack propagation angle defined on the basis of the Maximum Energy release rate criterion [30]. Material characteristics are reported in Table 2. In Fig. 13b, the initial mesh configuration used for the analyses is reported. According to the previous analyses, a relatively refined mesh is considered just around the crack tip region with characteristic length equal to ΔD =R ¼ 1=4 featuring Lagrange quadratic interpolation functions, whereas the remaining part of the structures presents a transition mesh with maximum characteristic length equal to ΔD =R ¼ 10=1. In this configuration, the total number of elements is equal to 1747 involving in 17684 DOFs. In order to achieve accuracy in the time integration a maximum time step equal to 5E-6 is adopted. In Fig. 16, the crack path

Fig. 11. RD law, loading configuation L3: crack tip speed as a function of the crack tip position (c/cr-XT/L) for the numerical simulations performed by using different mesh discretizations (M1, M2 and M3).

numerical formulations, such as SBFEM [28] and BEM method [29]. The geometric and loading configurations are reported in Fig. 13, in which a rectangular plate with length 140 mm, width 70 and thickness 15 mm is investigated. The loading scheme is based on symmetric compressive forces with respect to horizontal midspan axis acting along the end edges. In particular, the impact phenomena on the left and the right edges of the sample is simulated by considering uniformly distributed traction σ1 and σ2 , whose values are the ones collected during the ex­ periments (Fig. 14). This loading schematization refers to the mea­ surement system reported in Ref. [27], in which a striker bar produces impact phenomena by 2 bars located at the end edges of the specimen. An initial defect is assumed along the upper quarter length, with initial angle aligned with the horizontal axis. The present analysis can be classified as a “mixed phase simulation” case, according to the definition provided in Ref. [26], in which the time history of the crack tip speed is

Fig. 14. Time history of uniformly distributed traction on the end edges of the PMMA specimen with a hole.

Fig. 12. Synoptic representation of the mesh motion during the crack propagation for the M1 mesh configuration.

Fig. 13. PMMA specimen with a hole: geometry, loading and intial mesh configuation adopted for the simulations. 8

F. Fabbrocino et al.

Composites Part B 174 (2019) 107053

in agreement, since they denote the same trajectory of the crack tip front. From the analyses of the mesh configurations (Fig. 10), it is quite clear how the proposed model need to have accuracy only in a region close to the crack tip front, where a good prediction of the fracture variable is required. In the remaining region, a lower discretization is utilized, leading to a reduction of the computational costs of the nu­ merical model. However, as far the mesh motion introduces distortions in the mesh model due to crack evolution, in order to ensure accuracy in the FE model a remeshing procedure should be utilized. Although this may increase the computational cost of the model, the recourse to remeshing procedure is strongly reduced by the use of ALE approach, specifically by the regularization method. A synoptic representation of the deformed and remeshed configurations of the crack tip are reported in Fig. 17, whereas the crack tip position affected by remeshing pro­ cedure are presented in Fig. 16. Finally, comparisons in terms of dynamic SIFs time history deter­ mined between proposed model and numerical data taken from the literature are discussed (Fig. 18). In particular, mode I and mode II SIF are investigated. The evaluation of the SIFs is performed by the use of component separation method defined in Ref. [26]. The results show in both cases an unstable behavior, mostly for the mode II case, in which snap through phenomenon is observed. From the analysis of the results, instability phenomena are triggered by the initiation phase of the crack path. Subsequently, the SIFs for both mode I and mode II denote an oscillatory behavior. During the crack path the crack faces do not pro­ duce contact phenomena. Finally, a good agreement is observed be­ tween the predictions obtained by the proposed model and those taken

Fig. 15. Prescribed crack tip speed time history. Table 2 PMMA specimen with a hole: mechanical parameters. E ½MPa�

ν

μ ½MPa�

ρ ½kg m 3 �

cr ½m s

3300

0.42

1162

1180

1064

1

�

Fig. 16. Predicted crack path, comparisons with SBFEM [28] and BEM [29] models.

from the initial configuration is investigated. In particular, comparisons in terms of crack tip location during the crack growth between proposed model and numerical results, based SBFEM and BEM formulation, are presented. The results obtained by using different formulations are quite

Fig. 18. Dynamic SIFs time history: comparisons with SBFEM [28] and BEM [29] models.

Fig. 17. Synoptic representation of the mesh distorsion around the crack tip: (a) remeshed mesh i-th, (b) deformed mesh i-th, (c) remeshed mesh i-thþ1. 9

F. Fabbrocino et al.

Composites Part B 174 (2019) 107053

by existing numerical formulation available from the literature.

[9] Li F, Xingwen D. Mesh-dependence of material with softening behavior. Chin J Aeronaut 2010;23(1):46–53. [10] Murotani K, Yagawa G, Choi JB. Adaptive finite elements using hierarchical mesh and its application to crack propagation analysis. Comput Methods Appl Mech Eng 2013;253:1–14. [11] Bruno D, Greco F, Lonetti P. A fracture-ALE formulation to predict dynamic debonding in FRP strengthened concrete beams. Composites Part B: Engineering 2013;46:46–60. [12] Lee GH, Chung HJ, Choi CK. Adaptive crack propagation analysis with the elementfree Galerkin method. Int J Numer Methods Eng 2003;56(3):331–50. [13] Menouillard T, Belytschko T. Smoothed nodal forces for improved dynamic crack propagation modeling in XFEM. Int J Numer Methods Eng 2010;84(1):47–72. [14] Dirik H, Yalçinkaya T. Crack path and life prediction under mixed mode cyclic variable amplitude loading through XFEM. Int J Fatigue 2018;114:34–50. [15] Rodriguez-Ferran A, Perez-Foguet A, Huerta A. Arbitrary Lagrangian-Eulerian (ALE) formulation for hyperelastoplasticity. Int J Numer Methods Eng 2002;53(8): 1831–51. [16] Greco F, Leonetti L, Lonetti P. A novel approach based on ALE and delamination fracture mechanics for multilayered composite beams. Composites Part B: Engineering 2015;78:447–58. [17] Funari MF, Greco F, Lonetti P. Sandwich panels under interfacial debonding mechanisms. Compos Struct 2018;203:310–20. [18] Funari MF, Lonetti P. Initiation and evolution of debonding phenomena in layered structures. Theor Appl Fract Mech 2017;92:133–45. [19] Funari MF, Greco F, Lonetti P. Dynamic debonding in layered structures: a coupled ALE-cohesive approach. Frat Ed Integrit� a Strutt 2017;11(41):524–35. [20] Lonetti P. Dynamic propagation phenomena of multiple delaminations in composite structures. Comput Mater Sci 2010;48(3):563–75. [21] Bruno D, Greco F, Lonetti P. Dynamic mode I and mode II crack propagation in fiber reinforced composites. Mech Adv Mater Struct 2009;16(6):442–55. [22] Funari MF, Lonetti P, Spadea S. A crack growth strategy based on moving mesh method and fracture mechanics. Theor Appl Fract Mech 2019;102:103–15. [23] COMSOL. COMSOL Multiphysics® v. 5.2. Stockholm, Sweden: COMSOL AB; 2015. www.comsol.com. [24] Zhou F, Molinari JF, Shioya T. A rate-dependent cohesive model for simulating dynamic crack propagation in brittle materials. Eng Fract Mech 2005;72(9): 1383–410. [25] Nishioka T. Computational dynamic fracture mechanics. Int J Fract 1997;86(1–2): 127–59. [26] Nishioka T, Tokudome H, Kinoshita M. Dynamic fracture-path prediction in impact fracture phenomena using moving finite element method based on delaunay automatic mesh generation. Int J Solids Struct 2001;38(30–31):5273–301. [27] Gr� egoire D, et al. Dynamic crack propagation under mixed-mode loading comparison between experiments and X-FEM simulations. Int J Solids Struct 2007; 44(20):6517–34. [28] Ooi ET, et al. Dynamic crack propagation simulation with scaled boundary polygon elements and automatic remeshing technique. Eng Fract Mech 2013;106:1–21. [29] Fedelinski P. Computer modelling of dynamic fracture experiments. In: Key engineering materials; 2011. p. 113–25. [30] Funari MF, et al. A numerical model based on ALE formulation to predict crack propagation in sandwich structures. Frat Ed Integrit� a Strutt 2019;13(47):277–93.

5. Conclusions A numerical model based on fracture mechanics and moving mesh method is proposed. The analysis is presented in a dynamic framework, in which the crack tip path is defined in terms of proper crack function and kinking criterion. The numerical strategy is developed in the framework of ALE approach in which mesh regularization equations are introduced to simulate an arbitrary crack advance. The numerical modeling is validated by means of several comparisons under dynamic loading conditions. In particular, comparisons in terms of rate effects of the crack growth have shown the need to utilize rate dependent constitutive laws, in which the fracture toughness increases with the crack tip speed. The proposed modeling is validated by several com­ parisons with experimental and numerical data, which show accuracy and robustness of the numerical approach. Moreover, sensitivity ana­ lyses in terms of mesh dependence and time required for the solving procedure are also developed. Acknowledgements The work of F. Greco and P. Lonetti is supported by Italian Ministry of University and Research (P.R. I.N. National Grant 2017, Project code 2017J4EAYB; University of Calabria Research Unit). References [1] Ingraffea AR, Grigoriu M. Probabilistic fracture mechanics: a validation of predictive capability. CORNELL UNIV ITHACA NY DEPT OF STRUCTURAL ENGINEERING; 1990. [2] Ingraffea AR. Computational fracture mechanics. 2007. [3] Barbero E, Lonetti P. An inelastic damage model for fiber reinforced laminates. J Compos Mater 2002;36:941–62. [4] Camacho GT, Ortiz M. Computational modelling of impact damage in brittle materials. Int J Solids Struct 1996;33(20–22):2899–938. [5] Tawk I, et al. Composite delamination modelling using a multi-layered solid element. Compos Sci Technol 2010;70(2):207–14. [6] Unosson M, Olovsson L, Simonsson K. Failure modelling in finite element analyses: element erosion with crack-tip enhancement. Finite Elem Anal Des 2006;42(4): 283–97. [7] Rabinovitch O. Cohesive interface modeling of debonding failure in FRP strengthened beams. J. Eng. Mech. ASCE 2008;134(7):578–88. [8] de Borst R, Remmers JJC. Computational methods for debonding in composites. Mech. Response Compos. 2008;10:1–25.

10