Nội dung

Engineering Solid Mechanics 7 (2019) 263-278 Contents lists available at GrowingScience Engineering Solid Mechanics homepage: www.GrowingScience.com/esm Non-local averaging in composite micro-mechanical material models Sandeep Medikondaa and Ala Tabieia* a Department of Mechanical and Materials Engineering, University of Cincinnati, Cincinnati, OH, USA A R T I C L EI N F O Article history: Received 7 July, 2019 Accepted 6 August 2019 Available online 6 August 2019 Keywords: Unidirectional composites Micro-mechanical model Continuum damage mechanics Non-local damage Finite Element Method LS-DYNA ABSTRACT Strain-softening material models have conventionally had a pathological mesh sensitivity in finite element simulations and composite materials are indifferent to this problem. Spurious localization is inherent to the structural problem in strain-driven softening. This limitation is caused as the partial differential equations (that govern the structural problem) become ill-posed as the tangent modulus becomes negative, for which uniqueness of the solution with respect to the spatial discretization is lost. This causes the numerical results to unrealistically concentrate in a single layer of elements. A basic theory of overcoming mesh sensitivity is the non-local continuum theory. In this work, three different non-local models with isotropic weight function have been proposed and implemented to work in conjunction with a composite micro-mechanical material model. The effect of a weighing function in each of these formulations has been studied in detail. All three non-local formulations have been observed to produce a nice smeared effect of damage unlike the local damage models. © 2019 Growing Science Ltd. All rights reserved. 1. Introduction The exceptional ability of composite materials to absorb a high amount of kinetic energy during impact loadings, combined with a high strength-to-weight ratio, makes their use extremely rewarding in real world applications. This ability is due to the resistance offered to the growth of cracks which in turn results from the heterogeneous micro-structure of composite materials. The crack propagation is typically constrained at the boundaries between the constituent materials, which widely vary in their mechanical properties. However, propagation and coalescence of the cracks can eventually lead to a steady weakening of the material before a complete loss in the integrity of the structure. Continuum damage mechanics, allows us to describe the heterogeneous micro-processes, such as damage occurring because of the straining of materials as the macroscale behavior of the composite materials. The crack development in the material is typically quantified as a damage state depending on the failure mode, whose evolution causes a loss in the effective stress response or the stiffness of the structure. It is generally acknowledged that finite elements models with a fine mesh (smaller element size) yield more accurate results (Lee et al., 2011; More & Bindu, 2015). However, this is no longer the case for strain-softening materials. Numerical studies have shown that the results in strain-softening materials are * Corresponding author. E-mail addresses: tabieia@ucmail.uc.edu (A. Tabiei ) © 2019 Growing Science Ltd. All rights reserved. doi: 10.5267/j.esm.2019.8.002 264 essentially dependent on the finite element mesh (Sluys & De Borst, 1992; Sluys et al., 1993; de Borst & Verhoosel, 2018). This complication is well known (Belytschko et al., 1986; Lasry & Belytschko, 1988), and any type of strain localization phenomena (failure, in-elasticity, damage), will tend to localize in the smallest element of the finite element mesh. In other words, the smallest element in the mesh will tend to fail/erode before other elements. Also, finer the mesh, the energy dissipated by the numerical model decreases and tends to extremely low values, sometimes even to zero. Hence, the uniqueness of the solution with respect to the mesh size is lost, which is quite troubling from a numerical standpoint. This phenomenon can also be explained from a mathematical viewpoint. The partial differential equations of motion governing the structural problem will change characteristic type at the onset of softening, from hyperbolic to elliptic in dynamic problems and the opposite in static cases. This results in an ill-posed problem as the initial and boundary conditions for one class of equations are not applicable for the other. It should be noted that the strain localization isn’t purely a numerical phenomenon and it can be observed in certain experiments, such as the necking of a ductile material in tension. However, these strains are localized in some finite volume or a narrow band, the dimensions of which are characteristic of the material. Different remedies addressing this problem have been presented in the open literature and can be classified into the following categories/approaches: 1. Cohesive crack models/Cohesive zone models (CZM) (Barenblatt, 1962; Bazant, 2019; Dugdale, 1960): These models acknowledge the presence of a strong discontinuity and describe softening with the help of a traction-separation law (Park & Paulino, 2011). In these laws, as the surfaces separate traction increases until a maximum is reached and then subsequently reduces to zero indicating complete separation. The CZM model doesn’t necessarily represent any physical material, but describes the cohesive forces which occur when the material layers are being pulled apart. 2. Crack Band Model (Bažant & Oh, 1983; Pijaudier-Cabot et al., 1993): The idea behind this theory is to relate the specific volumetric energy, which is defined by the area underneath the stress-strain curve, with the fracture energy of the material. Hence, for an element that is representative of the characteristic process zone (and later a macroscopic crack) size, the product of the area under the stress-strain curve and the element size, corresponds to the fracture energy and is a material constant. Since, the process zone is represented by a band of elements with highly localized strain, the softening part of the stressstrain law is adjusted according to the element size. 3. Regularized Models: These are the generalized continuum theories which are applicable on the level of the constitutive equations. They typically use a characteristic length to prevent the localization of strain. The regularized models lead to a continuously differentiable displacement field and this results in a continuous strain field. Strain localization is typically manifested into a narrow band of high strain elements. These models are further divided into: a) Differential (gradient-enriched) formulations: This approach introduces higher special gradients into the equations of motion, making the stresses in the constitutive model dependent on strain gradients. The gradient formulations are further divided into Explicit (Aifantis, 1987) and Implicit gradient formulations (Peerlings et al., 1996). For a detailed account the reader is referred to the work of Simone (2007). At the core of these formulations, the damage is driven not only by an equivalent internal variable (typically strain) but also by its Laplacean (second order differential term resulting from Taylor’s series expansion). As long as the distribution of the internal variable is uniform, its Laplacean vanishes and the model responds similar to a local formulation. However, once the internal variable experiences its highest value, the curvature of its profile becomes negative and the presence of the Laplacean in the non-local equivalent internal variable results in a smaller value than the local one. The major difficulty with differential and especially the explicit gradient formulations is the difficulty to program them in the finite element method, mainly as the strain gradients are not normally computed or stored. S. Medikondaa and A. Tabiei / Engineering Solid Mechanics 7 (2019) 265 b) Integral formulations: This category of nonlocal models abandon the classical assumption of locality and define stress at a point to be dependent not only on the state variable (usually strain, damage) at that point but also on a distribution of the state variable in a vanishing region around the point (Bazant, 1986). In some of the early work, Pijaudier-Cabot and Bazant (1987) have shown that the use of integral non-local formulations as an efficient way of dealing with strain softening problems. In the following years, the use of nonlocal formulations has been extended to a wide range of models including softening plasticity, smeared crack models and microplane models. A detailed summary of these models can be found in the pioneering work of Bažant and Jirásek (2002). It must be noted that there are other remedies for addressing the mesh sensitivity problem of strainsoftening material models. One of them is the inclusion of rate-dependency (physical or numerical) in the evolution of damage or visco-plasticity in the plastic hardening/softening material models. This approach however, has been reported to be ineffective in explicit time integration scheme since the time step is extremely small to provide a stable solution (Lasry & Belytschko, 1988). Another approach is the use of generalized kinematic relations (dual equilibrium equations), which include either the continua with microstructure (e.g., Cosserat’s continuum mechanics) or the continua with nonlocal strain (e.g., nonlocal elasticity) (Lakes, 1995). Both these methods introduce a new level of complexity and sophistication by altering the equilibrium equations, which makes them difficult to implement in a FE program. The objective of the current work is to implement and study the effect of integral nonlocal damage formulations inside a three-dimensional strain-rate and pressure-dependent composite micro-mechanical material model. In the due process, the following contributions have been made to the current state of the art: 1. A computationally inexpensive non-local progressive damage strategy (using Weibull damage functions) has been proposed in an explicit framework that helps better control the distribution of damage and prevents strain localization typically seen in strain-softening composite materials. 2. Three modifications have been proposed to the implementation of the proposed non-local damage law along with a detailed study on the effect of a weighing function in each case. 2. Micro-mechanical Material Model The representative volume cell (RVC) used to develop the micro-mechanical relations is shown in Fig. 1. This RVC is the same as the one originally proposed by Pecknold and Rahman (1994) and further used in various micro-mechanical models by Tabiei and Aminjikarai (2009), Tabiei et al. (2005), Tabiei and Chen (2001) and Medikonda et al. (2017). However, for completeness the micro-mechanics relations are briefly discussed here. The unit cell is divided into three sub-cells: one fiber sub-cell, denoted as f, and two matrix sub-cells, denoted as MA and MB respectively. The effective stresses in the RVC are determined from the sub-cell values by combining 2 material parts: material part A consists of the fiber sub-cell f and the matrix sub-cell MA, and material part B consists of the remaining matrix MB using the iso-strain boundary conditions for all the directions. The dimensions of the unit cell are 1 × 1 unit square. The dimensions of the fiber and matrix sub-cells are denoted by Wf and Wm respectively as shown in Fig. 1. and defined as shown below: Wf  Vf ; Wm  1  W f where, Vf is the fiber volume fraction. (1) 266 Fig. 1. A representative volume cell of unidirectional fiber reinforced polymer composite Visco-plastic constitutive relations based on modified Bodner-Partom state variable model, initially proposed by Goldberg et al. (2005) and further enhanced by Zheng and Binienda (2008) have been used to represent the matrix sub-cells MA and MB. The full details of the model can be found in these references, however for completeness only the incremental form of those equations are given below.   1  Z 2 n   Sij  d   2 D0 exp        ij   dt    2   e    2 J 2   (2) I ij 2 I I deij deij ; 3 where, deeI  d ijI  d mI & d mI  (d11I  d 22I  d 33I ) / 3 (3) dZ  q( Z1  Z )deeI (4) d  q(1   )deeI (5) deeI  Stability against high strain increments has been ensured by implementing a 4-step Runge-Kutta integration scheme. The details of this algorithm can be found in Medikonda and Tabiei (2018). Fiber response is initially assumed to be elastic and transversely isotropic however the material is assumed to become orthotropic with damage evolution. As assumed by Tabiei and Babu (2009), the damages to the fibers are considered to be a result of only the direct loading in the current model as well.   f    C f   where, (6) f C  is the stiffness matrix which can be partitioned into direct and shear stress stiffness matrices f as follows: C   0C  0C  .   f fd 33 33 fs (7) S. Medikondaa and A. Tabiei / Engineering Solid Mechanics 7 (2019) 267 The direct stress stiffness matrix, should be symmetric and the following relationship should be obeyed:  ij  Ei  ji Ej (8) , i , j  1, 2,3 and i  j (no summation) The direct and shear stress stiffness matrices in terms of the properties of the fibers are: C   (1  d )E f (1  f ) 1 1 23  f f f ( 1    2    23 12 21 )     Symm.   C   G12f 0 0    f  G23 0   Symm. Gof12   fd fs (1  d1 )E1f (1  d 2 )E2f  12f  21f (1  23f  2 12f  21f ) (1  d 2 )E2f (1  12f  21f ) (1  23f )(1  23f  2 12f  21f ) (1  d1 )E1f (1  d 3 )E2f  12f  21f   (1  23f  2 12f  21f )  ( 23f  12f  21f )E2f (1  d 2 )(1  d 3 )   (1  23f )(1  23f  2 12f  21f )   (1  d 3 )E2f (1  12f  21f )  f f f f (1  23 )(1  23  2 12 21 )  (9) (10) It should be noted that the damage parameters d i , i  1, 2 , 3 , in the above relations follow progressive failure models and are discussed in the following section. The model can also account for the strain-rate dependency of certain glass fibers (Blazynski, 1987) using the following relations: G 12f  a G s s  G of12 (11) t  12 1 ss   log t 0  o (12)  dt   Once the stresses in all the constituent sub-cells have been obtained, they are then combined using the iso-strain boundary conditions to obtain the effective stresses of the RVC.  11RVC  W f2 11f  (1 W f2 ) 11R (13)  22RVC  Wf2 22f  (1 Wf2 ) 22R (14)  33RVC  Wf2 33f  (1  Wf2 ) 33R (15)  12RVC  W f Vs4 12f  (1  W f Vs4 )(1  d4 ) 12R (16)  23RVC  Wf Vs5 23f  (1  Wf Vs5 )(1  d5 ) 23R (17)  31RVC  Wf Vs 4 31f  (1  Wf Vs 4 )(1  d6 ) 31R (18) 268 Since, the use of iso-strain boundary conditions for shear isn’t quite realistic in a physical (Medikonda et al., 2017) sense, different ad-hoc shear volume fraction coefficients, Vs4 and Vs5, for the in-plane and transverse shear have been introduced and have values quite lower than the volume fraction of the fibers. Damage parameters d i , i  4 ,5,6 , represent the damages imposed on the matrix material and affect only the shear stresses of the resin. All damage functions and their subsequent non-local formulations will be discussed in the following sections. 3. Isotropic Nonlocal Formulation As previously mentioned in the introduction, the nonlocal approach consists of replacing a state variable by its nonlocal counterpart obtained by weighted averaging over a spatial neighborhood of each point under consideration. Hence in a domain field V, the corresponding nonlocal state variable is defined as: zn ( x )    ( x , ξ )z n (19) ( ξ ) dV ( ξ ) V where,  ( x , ξ ) is a given nonlocal operator. In an infinite body, the weight function depends only on the distance between the ‘source point’, ξ , and the ‘target’ point, x , and is given by the following relation:  ( x,ξ ) (20)  ( x,ξ )    ( x , ξ )dV (ξ ) V It should be noted that the weighing function  ( x , ξ ) is a monotonically decreasing non-negative function of the distance r  x  ξ . Typically, a Gaussian distribution is considered as the weight function. However, however, a wide range of weighing functions can be considered. The behavior of a few of these functions has been shown in Fig. 2. The Gauss distribution function is given by the following relation:  xξ  ( x , ξ )  exp    2 L2  2  ,   (21) where L, is a parameter reflecting the internal length of the nonlocal continuum and should be experimentally determined. A truncated quartic polynomial or the bell-shaped function is given by the following relation: 2 (22) 2 xξ .  (x,ξ )  1  L2 A typical green function used in the implicit gradient models, which has a weight function equivalent to the integral nonlocal models is given by the following relation:  ( x, ξ )   xξ 1 exp  1   2l L    .  (23) Lastly, a simple conical function and the 3-Parameter function based on the work of Pijaudier-Cabot and Bazant (1987) are given by the following relations:  ( x, ξ )  1  xξ L (24) S. Medikondaa and A. Tabiei / Engineering Solid Mechanics 7 (2019)   xξ  ( x , ξ )  1     L     p     269 (25) q The use of the length parameter L, has been considered to be ad-hoc by Schwer (2011), however most other authors have stated this to be a material parameter that has to be experimentally determined. Studying the effect of this parameter and the different types of weight functions discussed above, in the context of a micro-mechanical material model is a topic of interest and will be undertaken in the current work. For the results presented in Fig. 2, a value of 1 has been used except for the green’s function where a value of 0.5 has been used. In addition, a value of p  8 and q  2 has been used for the 3-parameter weight function. Fig. 2. Nonlocal Weight Functions 4. Numerical Implementation In the current section, we discuss the methodology used for the numerical implementation of a nonlocal model in a non-linear finite element code (LS-DYNA®) (Hallquist, 2006; Hallquist et al., 2014; Wang & Xia, 1998). In this context, Andrade et al. (2011) presented an algorithm for the numerical implementation of a Lemaitre-based ductile damage nonlocal model. A big drawback of this model is its demand for simultaneous access to all integration points in the finite element mesh, this would make the analysis painstakingly slow. However, building on the work of Tvergaard and Needleman (1995). and Cesar de Sa et al. (2010), Andrade et al. (2011) at a later stage proposed a nonlocal formulation that can be easily incorporated as a general feature for any user-defined material model (UMAT) in LS-DYNA. Hence, the drawback of accessing the neighboring integration points at once is overcome by adopting a strategy that saves and uses information of the state variable from a previous time step. The disadvantage of such as assumption is that it necessitates small time steps for enough accuracy. However, since the explicit time integration scheme of LS-DYNA naturally requires a very small-time step (less than the critical time step ( t  2 max ) in order to guarantee stable solutions this condition is easily met. It is 270 worth noting that LS-DYNA does offer the option of using nonlocal formulations through the keyword *MAT_NONLOCAL, which uses the 3-parameter weight function discussed in equation 25, however this option is limited to the use of very few elastoplastic models, nonetheless user-defined material models. Once the generic state variable that is desired to turned nonlocal has been decided upon (for e.g. z in * n1 , at eq. 19). Using the approximated nonlocal approach. The updated value of the nonlocal variable, z the current time is given as: zn*1  K nl zn1 (26) nl where, K is a nonlocal penalty factor defined as K nl  (27) zn , zn where, zn and zn are respectively the local and nonlocal values at the last converged time step. Adopting the Gaussian quadrature integration rule for Eq. (19), we get: (28) npgi zi   w j J j  ij zi , j 1 where, ij is the nonlocal operator that relates the Gauss points i and j located at global coordinates x and ξ respectively. In additions the quantities w j and J j are the Gaussian weights and Jacobian evaluated at Gauss point j . Lastly, npg i is the number of Gauss points that lie inside the nonlocal volume of interaction from point i . It must be noted that the factors w j , ij and J j are merely geometrical in nature and they depend on the finite element mesh itself rather that the constitutive model (UMAT). Hence, these factors only need to be calculated once, at the start of a simulation. The key part of the nl implementation lies in calculation of the nonlocal penalty factor K , which is then used to calculate the nonlocal value of the constitutive variable z , which in turn is selected depending on the choice of the model that is enhanced with nonlocality. Fig. 3. Schematic flowchart illustrating the implementation of nonlocal strategy in LS-DYNA (Andrade et al., 2011) S. Medikondaa and A. Tabiei / Engineering Solid Mechanics 7 (2019) 271 In the context of a micromechanical model discussed in the previous sections. Three slightly different formulations have been proposed, implemented and tested in the current work. 4.1 Formulation-1: The damage function is given by the following relation: d n 1 k   1  E damaged   k ij  max 1  exp    t|c  me   ij      mk    , d kn  ,    (29) where,  ij is the undamaged stress, ijt|c is the maximum strength and mk is the damage softening exponent in the corresponding directions (i.e., 11, 22, 33, 12, 23 and 31) and the sub-script k ranges from 1 to 6. The calculated damage from Eq. (29), is then used in reducing the stiffness of the lamina. Ekdamaged  (1  d kn 1 ) Ekun damaged (30) By using the non-local formulation described above, a non-local value of the damage variable can be calculated by using the damage from the previous time step ‘n’ in Eq. (28), this results in: npgi d kn   w j J j  ij d kn , (31) j 1 which can then be used in calculating the nonlocal penalty factor and the updated nonlocal damage variable using the following equations: d kn K  n dk (32) d k( n 1)*  K nl d kn1 (33) nl Lastly, instead of the local value of the damage the updated nonlocal value can be used in reducing the stiffness of lamina. Ekdamaged  (1  d k( n 1)* ) Ekun damaged (34) The reduced values of the stiffness are then used to calculate the effective stresses of the representative volume cell discussed in the previous chapters. 4.2 Formulation-2: In the current formulation strain is non-localized instead of damage which leads to the following nonlocal strain calculation. npgi  ijn   w j J j  ij ijn (35) j 1 Based on which the nonlocal penalty factor is calculated using the following relation: K nl   ijn  ijn (36) 272 and subsequently this factor is used in the damage function given by the following relation: d n 1 k   1  E damaged K nl   k ij  max 1  exp    t|c k  me       mk    , d kn     (37) 4.3 Formulation-3: Like Formulation-2, the non-local strain is calculated based on equation (35). However, instead of calculating the penalty factor this calculated non-local strain is directly used in the damage formulation as shown below: d n 1 k   1  E damaged   k ij  max 1  exp    t|c me   k       mk     , d kn     (38) 5. Results and Discussion To test the non-local technique, it has been implemented along with the micro-mechanical material model discussed in the previous chapters and applied on a tensile dog-bone specimen, commonly used in the experimental determination of the properties of composites. The dimensions of the specimen are shown in Fig. 4. Fig. 4. Tensile Dog-bone Specimen (Andrade et al., 2011) Fig. 5 shows the 3 meshes that have been built for the specimen with different mesh sizes. It should be noted that the mesh size has been only significant increased for the curved part of the specimen since this is the primary area of interest. The symmetry in the specimen has been taken into consideration and hence only a quarter of the model has been modeled to reduce computational effort. Additionally, it must be noted that only one layer of elements with a 0-degree fiber orientation representing a single ply have been numerically modelled. Fig. 5. Tensile Dog-bone Specimen with different Fig. 6. Local Damage in the Tensile Dog-bone mesh sizes Specimen for different mesh sizes