Multiscale progressive damage analysis of CFRP composites using a mechanics based constitutive relation

Composite Structures 235 (2020) 111759

Contents lists available at ScienceDirect

Composite Structures journal homepage: www.elsevier.com/locate/compstruct

Multiscale progressive damage analysis of CFRP composites using a mechanics based constitutive relation

T

Rudraprasad Bhattacharyya, Prodyot K. Basu

⁎

Vanderbilt University, VU Station B#351831, 2301 Vanderbilt Place, Nashville, TN 37235, USA

ARTICLE INFO

ABSTRACT

Keywords: Progressive damage analysis Composite laminate Multiscale modeling Damage mechanics

A continuum damage mechanics-based material model for the constituent phases of a carbon fiber reinforced polymer composite material is presented. This constitutive model is applied to Eigenstrain based Reduced-order Homogenization framework for multiscale modeling. A key contribution is that all model parameters are directly or indirectly calibrated using experimental data and these model parameters can be readily updated by describing the evolution of damage as a function of the strain state. The constitutive model for polymer matrix captures mixed mode from pure uniaxial to pure shear loading without using any conventional failure criterion. The accuracy of the constitutive model is demonstrated by comparing experimental results of unnotched and open-hole laminates. In addition, a demonstration is presented to alleviate spurious residual stiffness appearing from the reduction of model order of the microstructure.

1. Introduction The design and certification of aerospace composite structures require the extensive use of experimental data for a range of specimens, from coupons to full-scale prototypes. Apart from experimental testing being expensive and time-consuming, this approach is prone to overconservative design. In the case of composites, chosen for its superior macroscopic properties, the lack of a full understanding of damage initiation and growth at the microscopic level culminating into failure at the macroscopic or structural level precludes the use of the damage tolerant approaches commonly used in the design of metallic structures. As a result, conventional composite airframe designs cannot fully capitalize on the weight and performance gains that are possible with composites. As the potential benefits of incorporating predictive failure modeling and simulation methodologies in aircraft design, certification, and maintenance processes are being recognized, there is significant effort aimed at improving predictive computational methods, as exemplified by the ongoing efforts at U.S. DoD [1] and NASA [2,3]. A properly validated predictive modeling scheme is critical to rapid and economical aircraft design and certification processes with new composite materials, including novel nano-engineered composites. Progressive damage analysis (PDA) is a broad term applied to modeling methodologies that allow for the prediction of the initiation and evolution of damage. Several material constitutive models have been proposed in the existing literature, namely, Pinho et al. [4,5],

⁎

Ridha et al. [6], Leone et al. [7], Su et al. [8]. In addition to continuum damage models, computational schemes based on cohesive zone approach are also presented in the literature, such as Higuchi et al. [9], van Dongen et al. [10]. Progressive Damage Analysis in the context of multiscale modeling is presented in Laurin et al. [11], Pineda et al. [12], and Massarwa et al. [13]. The proposed Eigenstrain-based Reduced-order Homogenization Method (ERHM) modeling approach presented herein is in line with the hierarchical multiscale method described by Fish [14]. In this method, the nucleation and propagation of damage within a composite are tracked at the level of constituent material using continuum damage model, without applying any phenomenological failure criterion at the macroscale. In ERHM, the microstructural information like localization operators for mapping micro-scale behavior onto macro-scale properties, concentration tensors used in the homogenization process, and influence functions characterizing the influence of macroscopic displacements on microscopic response, are precomputed before undertaking the progressive damage analysis at the macroscale. The microstructural analyses are performed over a Representative Volume Element (RVE), which are coupled with structural analysis at the continuum scale. Computational efficiency of homogenization is reduced by developing a reduced-order microstructural geometry. The microstructural analyses are concurrently performed during progressive damage analysis. ERHM can account for progressive debonding between fiber and matrix at the scale of the microstructure and is equipped with an adaptive model

Corresponding author. E-mail address: [email protected] (P.K. Basu).

https://doi.org/10.1016/j.compstruct.2019.111759 Received 3 July 2019; Received in revised form 25 November 2019; Accepted 28 November 2019 Available online 04 December 2019 0263-8223/ © 2019 Elsevier Ltd. All rights reserved.

Composite Structures 235 (2020) 111759

R. Bhattacharyya and P.K. Basu

Fig. 1. Schematic representation of Macro- and Micro-scale problem in CFRP.

improvement capability to hierarchically increase model fidelity, such as those used by Oskay et al. [15] and Crouch et al. [16]. Crouch et al. [17] and Bogdanor et al. [18] used a symmetric version of ERHM (SERHM) to predict the failure of open-hole specimens under static loading. In all these studies, an arctangent type damage potential function was used. The major drawback of using this function entails lesser control in the softening regime due to the asymptotic nature of the potential near the failure strain. This leads to artificially enhanced ductility in predictions at the continuum scale. Moreover, the calibration of damage model parameters is more qualitative than mechanics based. Composite materials under shear loading exhibit a nonlinear stressstrain response prior to failure. Several evidences of shear nonlinearity in laminated composites embodying with high strength fibers embedded in epoxy resin can be found in Soutis et al. [19], Gilat et al. [20], Littell et al. [21], Park et al. [22], and Salavatian et al. [23]. This nonlinear plane shear response of the composite is widely known, and there are several approaches to consider this effect in material modeling, e.g., those by Paepegem et al. [24] and Lomakin et al. [25]. Computational perspectives are addressed in Paepegem et al. [26] and Fedulov et al. [27]. In contrast, brittle failure is observed in composites under pure uniaxial loading. In damage mechanics based models, the ductility of a material is controlled by model parameters of damage evolution. A novel weighting approach for the damage evolution model is proposed in this paper, to account for the differences in damage evolution of composites under purely normal and purely shear loads. Earlier, in Bogdanor et al. [18], the range of strain triaxiality was varied from pure hydrostatic to pure shear loading. No experimental evidence is available to validate the general shape of the failure surface employed. This can be attributed to the fact that the constitutive model represents the behavior of in-situ resin. To address the stated concerns, a new damage model is proposed. This model is incorporated into the multiscale modeling framework, verified and validated using numerical simulation of composite laminates in different configurations. This damage model switches from an arctangent expression into a simplified linear function in the post-ultimate or softening stress state. The model is designed to ensure satisfaction of the following criteria.

presented in Section 2. An elastic-damage constitutive behavior for the composite constituents is described in Section 3. The implementation of the damage model in the multiscale modeling framework is considered in Section 4. Incorporation of impotent eigen-modes within the multiscale framework is shown in Section 5. Calibration and validation of the damage model using numerical results are addressed in Section 6. 2. Multiscale modeling framework Consider a three-dimensional macroscopic fiber composite domain, 3 , made of a heterogeneous, periodically repeating representative 3 . The representative volume consists of two constituent volume, material phases, namely, carbon fiber and polymer matrix. The multiscale structure of the material is schematically illustrated in Fig. 1. Due to heterogeneity in fiber composites, the displacement fields fluctuate at the scale of the microstructure, . A scaled coordinate system, y = x/ , is introduced to capture the scale of material hetero1. The geneity. The scaling parameter, , satisfies the condition 0 < representative volume, , is parameterized in terms of y coordinate. The macroscale stress ( ) and strain ( ) are obtained by spatial averaging of the stresses and strains, as shown in Eqs. (1) and (2).

x, t =

x, t =

1

1

x, y, t d y

(1)

x, y, t d y

(2)

in which is the volume of the representative volume element. The equilibrium equations at two different scales, namely, O ( 1) and O (1) are expressed as:

O(

1):

y ·[C (x ,

y): ( ¯ (x, t ) +

1 y u (x ,

y, t )

µ (x, y, t ))] = 0 (3)

O (1):

x · (x , y , t ) = 0

(4)

in which the tensor of elastic moduli, C , is considered to vary as a function of the microscopic coordinate (y) only, µ are the inelastic strains in the microstructure due to damage and u1 is the locally fluctuating displacement. Eq. (3) constitutes the microscale equilibrium equation applied over the domain of the representative volume, . Eq. (4) is the result of averaging over to obtain the macroscopic equilibrium equation defined over the problem domain, .

• Ability to predict response under the conditions of pure shear, pure uniaxial (tension or compression) and mixed-mode loading. • Ability to smoothly capture the triaxial state of stress. • Ability to consistently control the softening regime and material failure.

2.1. Eigenstrain based reduced order modeling

The mathematical form of the damage model is verified for both isotropic and transversely isotropic damage susceptible elastic constituent materials. In addition to the development of a new constitutive model, the multiscale framework is improved upon to account for post failure stiffness and strength losses. A brief description of the multiscale modeling framework is

Implementation of mathematical homogenization using Eqs. (3) and (4) leads to the standard computational homogenization (also known as FE2) scheme proposed by Feyel et al. [28]. As the direct implementation of FE2 approach is computationally expensive, an eigenstrain based reduced order model is developed herein following Oskay et al. [15] and Crouch et al. [16]. For completeness, a concise overview of the 2

Composite Structures 235 (2020) 111759

R. Bhattacharyya and P.K. Basu

reduced-order modeling approach for computational homogenization is put forward. In a heterogeneous microstructure, the microscale displacement field () is defined in terms of the macroscopic applied strains, ( ) and the inelastic strains at the microscale, ( µ ) as below.

u1 x, y , t = H (y):

x, t +

h y, y : µ x , y , t d y ;

~=

n

s yu

( )

P (y): µ (

=

)

x, t

=1

where,

y, y

( )

P (y) =

g y , y N ( ) (y ) d y ;

( )

y, y

(14)

(5) where, H and h are third order tensors representing the elastic influence function and damage induced influence function, respectively. These influence functions are based on the elastic behavior of the microstructure, both in the absence and presence of damage. For reduced-order modeling, the microstructure response is based on a reduced number of degrees of freedom. By precomputing the coefficient tensors for displacement in the representative volume, the reduced-order modeling approach economizes the computational cost, as per Oskay et al. [15]. The representative volume domain, , is partitioned into n non-overlapping subdomains, ( ) , referred to as failing part or damaged part, with = 1, 2, …n indicating the part number. Mathematically, it can be expressed as below. n

( );

( )

( )

2.2. Macroscale boundary value problem The macroscale equilibrium and kinematic equations as well as boundary conditions under quasi-static conditions are as follows.

· (x, t ) + b (x) = 0; t

(x , t ) =

( )

n

( )

x, t

n

N ( ) (y) µ (

)

x, t

(8)

n

=

B(

)

F(

=0

P( ) : µ(

+

= 1, 2, …, n

)

)

P(

(11)

)

),

F(

),

L( )

( )

( )

P (y): C (y): A (y) N ( ) d y ( )

=

=

(19)

=1

=

L( ) =

( )

( )

P (y): C (y): P N ( ) d y

1 ( )

1

C (y): A (y) N ( ) d y

P( )

are coefficient tensors (20) (21) (22)

( )

( )

C (y): P (y) N ( ) d y

(23)

Since the microstructural morphology and elastic constants of the constituents are treated as deterministic, the coefficient tensors are also deterministic. The homogenized macroscopic stress in Eq. (19) constitutes stress update in the solution for the macroscale boundary value problem.

(12)

where I is the fourth order identity tensor and G is the elastic polarization tensor, evaluated as: s yH

)

In Eqs. (18) and (19), and due to Crouch et al. [16] defined as below.

(9)

The fourth-order elastic strain concentration tensor, A , is expressed as:

G=

L( ): B(

where, the fourth-order part damage polarization tensor, g , is given by

A = I + G (y)

µ(

n ( ))

(1

(10)

y)

):

=1

Here, n is the number of damaged parts in . The homogenized macroscopic stress is computed as a function of part averaged damage and eigenstrain coefficients, as shown below.

(x, y , t )) g y, y C (y)·[A (y): (x, t ) + ~ (x, y, t )] d y = 0

s y h (y ,

F(

+

(18)

In Eq. (5), the microscale displacement field is decomposed into linear and damage induced components as in Crouch et al. [16]. Consequently, the governing equation of the representative volume takes the form:

g (y , y ) =

):

=1

=1

(1

B(

=1

n

µ x , y, t =

( ))

(1

n

N ( ) (y)

(17)

t

Given the macroscale strain state, ¯ , and damage measure, ( ) , in each damaged part, , of the microstructure defined by the evolution equations for part-averaged damage, the evolution of microscopic damage is described by continuum damage mechanics based on the history dependent strain state in each damaged part and is stored as an internal state variable in the macroscale model. The unknown eigenstrains in each damaged part, µ ( ) , are obtained using Eq. (18).

In this equation, the subscripts and identify the constituent phases in a subdomain ( ) or a damaged part, . Implementation of the reduced order model for a unidirectional fiber composite is described in Section 4. Shape functions, N ( ) , are defined as piecewise constant functions, which form a partition of unity in the characteristic volume (i.e. ( ) and 0 otherwise). The reduced basis representation N ( ) = 1 if y for damage fields, and damage induced inelastic fields, µ , of nth order reduced basis are represented as below.

=1

· n = t (x , t ) x

u

2.3. Microscale boundary value problem

(7)

x , y, t =

(16)

where, is the homogenized Cauchy stress field, b is the body force. The spatial and temporal coordinates are x and t. u¯ is the macroscale displacement field. The displacement and traction boundaries are denoted by u and t , respectively. ·(·) and s (·) are the divergence and symmetric gradient operators, respectively.

The reduced-order model is partitioned to avoid parts with no intersecting constituent phase, as shown in Eq. (7). The parts of the representative volume are partitioned to enable grouping of the similar regions of the microstructure corresponding to the macroscopic failure modes. ( )

(15)

t)

u¯ = u (x, t ) x

(6)

=1

su ¯ (x ,

[0, t f ]

3. The constitutive model

(13)

The part damage induced strain, ~ , can be explicitly expressed as:

The evolution of the reduced-order representative volume problem 3

Composite Structures 235 (2020) 111759

R. Bhattacharyya and P.K. Basu

in S-ERHM is governed by the constitutive relationship between the average microstructure displacement variables. For this purpose, a rateindependent continuum damage evolution model is proposed. The damage evolution expressions are used to compute the failure strain within the context of the damage model employed in the multiscale formulation described in Section 2. In this paper, the constituent material properties of carbon fiber and polymer matrix are assumed as isotropic and transverse isotropic. Employing the conventional formulation for continuum damage mechanics, the stress-strain relationship in each reduced-order model damage part is expressed as: ( )

( )

(x , t ) = (1

(x, t )) C ( ) :

( ) (x ,

t)

evolution functions defined by Eq. (29):

0

(

=

1 ( (F 2

( ) T

)

) : C( ): F

( )

( )

F

x, t =

h( ) =

1 c(

if )

)

0

0

h 2(

0

0

( )

(25)

( )

>0

otherwise

)

kbold =

=

( )

(t ) ;

(

( )

(t ))

0

( )

(1 + A f( ) )] + (1 + A f( ) ) if

( ) tr

<

( )

2

( ) tr )

(

if

( )

( ) tr

<

( ) f

( ) f

(

max max

2

+

max

0, 1

)

(30)

where, max is the maximum engineering shear strain and max is the maximum absolute principal strain. For the sake of simplicity, the superscript ( ) has been dropped. The parameters controlling damage evolution in any damaged part are weighted as a function of maximum principal strain ( max ) and maximum shear strain ( max ) , as defined below.

(27)

in which, c ( ) is a material parameter that represents the damage anisotropy of the tensile and compressive loading in principal directions. ( )

<

One of the primary objectives of the constitutive model for the polymer matrix is to capture the strain triaxiality effect present under pure shear loading. In this context, the concept of triaxiality factor, kbold , was used by Bogdanor et al. [18] in the form defined in Eq. (30).

(26)

for = 1, 2, 3

( ) 0

3.1. Modeling of hardening regime

0 h3(

if

+ atan( ( ) )

in which, the damage evolution parameters, ( ) and ( ) , control the nonlinearity and strength of a damage part , depending upon the constituent material. 0( ) is the value of the damage equivalent strain at the initiation threshold for damage, tr( ) is the damage equivalent strain value at the transition point, and f( ) is the damage equivalent strain value at failure. So, in addition to ( ) and ( ) , this proposed model is characterized by these three damage equivalent strain parameters. A f( ) is the slope of the softening regime.

0 )

=

( ) + atan( ( ) )

( ) 0

( )

(29)

(24)

( )

h1(

( ) 0

( )

1

where, F and represent the strain weighting matrix and the principal strain vector within the damaged part . The superscript T indicates the transpose operation. The strain weighting matrix is computed as: ( )

( ))

( ) tr [ ( )

in which, ( ) and ( ) denote the part-averaged stress and strain tensors for each damage part, respectively; and C( ) is the tensor of elastic moduli for the constituent (i.e., fiber or matrix) that occupies the damaged part, . Damage accumulation in each constituent damaged part is driven by the phase damage equivalent strain, ( ) . ( )

( )

atan

if

(28)

max

= max{

In the present study, a new damage potential is proposed. The hardening segment of the constitutive model is defined by an arctangent evolution function, whereas the softening segment is defined to follow a linear law. This damage model is referred to as Arctangent Hardening Linear Softening (AHLS) model. A schematic stress-strain diagram is presented in Fig. 2. The transition strain, (tr ) , for the damaged part is defined as a strain corresponding to the ultimate stress or the strain in close proximity of the strain corresponding to the ultimate stress. This parameter is derived from transition damage equivalent strain tr( ) , which is an algorithmic parameter. The new damage potential function is described by two different

max

=

2

1

,

2

max{ }

,

(31)

}

3

min{ }

(32)

2

The numerical value of is 0 for pure hydrostatic loading, and for pure shear loading the value is 1. Since the experimental data used for damage model calibration involve uniaxial as well as shear loading on different laminate configurations, the calibration of damage model parameters using Eq. (30) is inconvenient. To overcome this problem, the following alternative expression for kb is proposed.

kbold

kb =

(3 + ) (1 )

(

max max

2

+

max

)

2(1 + ) (3 + ) +

(33)

Here, is the Poisson’s ratio of the polymer matrix, and · + denotes the Macaulay bracket. The numerical value of kb for pure uniaxial loading is 0 and for pure shear loading it is 1. This new expression for triaxiality factor is expected to ease parametric calibration based on experimental data. The derivation of this new triaxiality factor is shown in A. The variables of damage evolution under an arbitrary loading state are expressed as a function of this triaxiality factor. The parameters are calibrated based on pure shear (subscript S) and pure uniaxial (subscript N) load states, as shown in Eqs. (39)–(41). Based on the triaxiality state of the stress, the weighting factor kb is evaluated using Eq. (33). Then the damage evolution parameters, and ( ) , are weighted using the triaxiality factor, kb . The weighting parameter is specific to a given material. These two triaxiality parameters dictate the values of parameters, and , which together control the extent of ductility in the hardening regime as well as the ultimate strength of the material.

Fig. 2. A typical representation of the stress-strain diagram for constituent. 4

Composite Structures 235 (2020) 111759

R. Bhattacharyya and P.K. Basu

= (kb)

S

+ (1

kb )

N

(34)

= (k b )

S

+ (1

kb )

N

(35)

Young’s modulus, 12 is the longitudinal Poisson’s ratio and transverse Poisson’s ratio.

=

(T) (·)

E 2(1 + )(1

2 )

(S) (·)

=

(S) (·)

E 2(1 + )(1

2 )

(C) (·)

=

(C) (·)

E 2(1 + )(1

2 )

1

+ 2 2 (c 2

(1 + c 2 )(1

c 2 (1

)

2c ) =

(T) (·)

(m ) T

2 c =

(S) (·)

(m ) S

) + 2 2 (1

2c ) =

(C) (·)

The underlying assumption in the damage potential model of Eq. (29) is that the softening behavior of each reduced-order part can be modeled using a linear stress-strain law, which allows nonlinear damage evolution within the hardening regime. The damage evolution law derived here corresponds to linear softening and follows C 0 continuity. The spatial (not temporal) stress rate tensor in terms of strainrate tensor is computed using Eq. (48).

= (1

The solutions of the differential equations for a polymer matrix constituent under pure uniaxial and pure shear loading are presented in Eqs. (49) and (50), respectively.

(m ) C

(38) In Eqs. (36)–(38), E is the Young’s modulus of the constituent polymer matrix material. The subscript (·) is a generalized representation of the damage initiation, transition and failure of the constituent matrix. The damage equivalent strain parameters corresponding to damage initiation, transition and failure under pure uniaxial loading are denoted by (N) (N) and f(N) , respectively. Similarly, damage equivalent strain 0 , tr parameters corresponding to damage initiation, transition and failure under pure shear loading are denoted by 0(S) , tr(S) and f(S) , respectively. Under a mixed mode state of strain loading, the parameters for damage initiation, transition and failure are evaluated using Eqs. (39)–(41). (S) 0

+ (1

(k b ) )

(N) 0

(39)

tr

= (k b )

(S) tr

+ (1

(kb) )

(N) tr

(40)

f

= (k b )

(S) f

+ (1

(k b ) )

(N) f

(41)

(C) (·)

=

(T) (·)

=

(C) (·)

C33 + 2c 2

2 31 (C11

+ C12 ) 2

4c

31 C13

c 2C33 + 2

2 31 (C11

+ C12) 2

4c

31 C13

=

(T) (·)

(f ) T

=

(C) (·)

(f ) C

C12 =

C13 = C33 =

1

2 12

1

2 12

1

2 12

1

2 12

E11 (1 2 13

13 31) 31

2

12 13 31

E11 ( 12 + 2 13 31

13 31)

E11 ( 31 + 2 13 31

12 31)

E33 (1 2 13

2

2

12 13 31

12 13 31

2 12 ) 31

2

12 13 31

=

(S)

=

A f(N) =

(N) tr

(S) tr

[ (

[ (

(1 + Af(N) )] + (1 + Af(N) )

(N) tr )

(S) tr )

(1 + Af(S) )] + (1 + Af(S) )

(49) (50)

max

E(

(N) f

(N) tr )

;

Af(S) =

max (1 (S) f

E(

+ ) (S) tr )

where, max and max are ultimate stress under pure uniaxial and pure shear loading. The elastic tensor in the softening regime is based on the weighted sum of the elastic contributions in pure uniaxial and pure shear loading, as shown in Eq. (51):

Af = kb A f(S) + (1

kb) Af(N)

(51)

The constituent carbon fiber is assumed to demonstrate brittle failure under uniaxial loading. The mathematical expression for damage evolution in fiber is same as Eq. (49), with

A f(N) =

(42)

max

C33 (

f

tr )

(52)

4. Numerical implementation Progressive damage evaluation of composites involve analyses at two different scales – micro and macro. These two analyses are undertaken separately. The Graphical User Interface (GUI) in the commercially available software Abaqus 6.14 is used for creating the morphology of the representative volume and the macroscale coupons. The microscale analysis is based on the in-house software CoefTensCompute (CTC), originally developed by Crouch et al. [16]. The macroscale analysis was undertaken using the user-defined material model subroutine (UMAT) capability of Abaqus. Details of microscale and macroscale analyses are given in Sections 4.1, 4.2.

(43)

where, C11, C12, C13 and C33 are components of the tensor of elastic moduli for transversely isotropic carbon fiber. These components are given by:

C11 =

(N)

In Eqs. (49) and (50),

By extending the AHLS damage model to the case of transversely isotropic fiber, the expression for damage equivalent strain parameters is given by Eqs. (42) and (43). (T) (·)

(48)

C:

1. Damage parameter reaches unity, (i.e. complete loss of load carrying capacity) at failure strain ( f ) = 1. 2. Stress remains continuous at the transition point: ( tr ) = ( tr+)

(37)

= (k b )

) C:

In order to ensure continuity, the damage model must satisfy the following constraints:

(36)

0

is the

3.2. Modeling of softening regime

Here, N, N are the damage model parameters for pure uniaxial loading. On the other hand, S, S denote the damage model parameters for pure shear loading. Note that in these equations, the subscript N becomes T for tension and C for compression. For a specific strain state, the direction of axial loading is evaluated first and the parameters are selected accordingly. These four damage model parameters can be calibrated directly with the experiment data. In this proposed damage model, the damage evolution for each loading condition (either, pure uniaxial or, pure shear) depends on the corresponding three damage equivalent strain parameters. These parameters can be obtained using Eqs. (36) and (37). The details of formulation and evaluation of the damage equivalent strain parameters are given in B. (T) (·)

31

(44) (45)

4.1. Microscale analysis

(46)

The microscale analysis of a composite representative volume is performed prior to the simulation as a preprocessing step. The geometry of the microstructure and corresponding reduced-order model parts are generated first.

(47)

Here, E33 is the longitudinal Young’s modulus, E11 is the transverse 5

Composite Structures 235 (2020) 111759

R. Bhattacharyya and P.K. Basu

Fig. 3. Partitioning of the representative volume and Reduced Order Model.

Fig. 4. Pictorial representation of the concept of periodicity.

4.1.1. Creating reduced order model morphology Inputs to microstructural analysis are morphology of the representative volume (e.g. cell type, fiber volume fraction) and elastic parameters of constituent materials. In this study, a cubic packed fiber composite microstructure is assumed. Accordingly, a representative cubic volume consisting of the matrix and fiber material phases is modeled to represent the microstructure of the unidirectional fiber reinforced composite with approximately 65% fiber volume fraction. The microstructure domain of size 10 µ m× 10 µ m× 10 µm is discretized with 1640 hexahedron (also referred to as C3D8) elements and 2035 nodes. This representative volume is then partitioned, as shown in Fig. 3, into 4 distinct damaged parts corresponding to the macroscale failure modes of the composite, namely, fiber failure, transverse matrix cracking, and delamination. The ‘interaction damage’ part represents the maximum damage accumulated by either, transverse matrix cracking or, delamination. The ‘fiber damage’ part consists of 960 elements and 1243 nodes. The ‘matrix cracking’ and ‘delamination’ damage parts consist of 320 elements and 594 nodes. The ‘interaction part’ part consists of 40 elements and 176 nodes. Unlike Crouch et al. [16] and Bogdanor et al. [18], perfect bonding is assumed at the interface and no additional interfacial properties between the fiber and matrix within the representative volume are considered. For microscale analysis, all eight vertex nodes are assumed as fixed in all three directions. Then the periodicities in the representative volume are identified. The periodicity of displacements on the boundaries

of the representative volume is allowed for by identifying edge periodicity and face periodicity. For the sake of computational convenience, the two periodicity conditions are enforced as below. 1. Edge periodicity: In the representative volume element there are 4 edges parallel to each spatial coordinate. Considering one spatial direction (say, Z-direction), one among the four edges is designated as the ‘master edge’ and the remaining three edges in the same direction are designated as ‘slave edges’. An injective mapping is performed from ‘master edge nodes’ to ‘slave edge nodes’. In Fig. 4a, the nine nodes marked with solid blue circles on the edge highlighted with a blue line which is parallel to the Z-axis are termed as master nodes. The remaining three edges parallel to z-axis marked with nine open red circles are termed as slave edge nodes. Likewise, master and slave nodes are identified in the remaining eight edges, parallel to x- and y-axes. The edges parallel to x-axis is marked in red and the ones parallel to y-axis in pink. Master-slave node pairing is also done for the edges parallel to X and Y directions. 2. Face periodicity: In face periodicity, 3 pairs of faces parallel to X-Y, Y-Z and Z-X planes are considered. For example, for Z- plane (or, X-Y plane) the face periodicity is shown in Fig. 4b. All the nodes on the back faces of the cube are considered as ‘face master nodes’. Corresponding slave nodes on the front face in XY plane are marked with orange circles. Similar face periodicities are identified for the faces corresponding to Y- and Z- planes. In mathematical sense, this 6

Composite Structures 235 (2020) 111759

R. Bhattacharyya and P.K. Basu

again is injective mapping.

III) The solution of the micro-problem is performed at each integration point of the macroscale finite element. The precomputed coefficient tensors expedited the macrostructural analysis, as on-thefly evaluation of the microscale model could be undertaken for the reduced order model (with reduced number of parts), without the need for using fully resolved microstructure. IV) At each incremental step, the triaxiality parameter, kb , is evaluated for all three reduced order damaged parts, namely ‘transverse matrix cracking’, ‘delamination’, and ‘fiber cracking’. Accordingly, appropriate damage model parameter is evaluated to capture the damage state, using Eqs. (39)–(41). V) The shape parameters are weighted using Eqs. (34) and (35). Once the material starts to soften, damage evolution is governed by the parameters, tr and Af , only, as with Eq. (29). VI) Since damage is considered directly at the scale of the constituents, damaged parts are a consequence of the microscale response, which is explicitly evaluated within the multiscale model. VII) Ensuring convergence in the solution for the reduced order model problem for each macroscale finite element, the UMAT computes the homogenized secant stiffness tensor and stress at each integration point based on the macroscale strain given by the finite element solver. VIII) At each incremental step, for a macroscale element the damages and damage-induced strains within each reduced order part are passed on to Abaqus as state variables.

These definitions ensured the periodicity of the representative volume during the computation of coefficient tensor. In this step, Abaqus compatible input file (Input_filename.inp) is converted into a ‘CTC’ compatible file (MicroInput.inp) using the Python script (named as AbaInp2MicroInp_conversion.py). The node and element sets are predefined in the Abaqus compatible input file. All the node and element descriptions are processed first. Next, the assembled element and node sets are identified. Finally, the representative volume is discretized with 1640 finite elements, and four damage parts, as shown in Fig. 3, and efficiently analyzed as a reduced order model. 4.1.2. Computation of coefficient tensors for reduced order model The scheme for the computation of coefficient tensors of the composite reduced order model using CTC is described as follows. I) The geometrical features of the representative volume domain and elastic material parameters are input. II) Elastic and part damage polarization functions are computed using Eqs. (11) and (13). Coefficient tensors A and P are also computed in this step using Eqs. (12) and (14). III) Partitioning of damaged part for the representative volume is constructed. IV) Coefficient tensors B, F, L and P are obtained using Eqs. (20)–(23). V) Finally, all the coefficient tensors are stored in an output file.

Python script is used to post-process the information in the Abaqus output data files generated from the numerical simulation. Information on stress and strain are extracted from this file to create the stress-strain plots. In this study, six noded tetrahedron (C3D6) and eight noded hexahedron (C3D8R) elements with reduced integration points are used for domain discretization. In the case of angle plies inside a laminate, the mesh is aligned in the orientation of the fiber to avoid mesh bias in the solution. In Fig. 6, the computational scheme followed in Section 4.1, 4.2 is shown as a flowchart.

4.2. Macroscale analysis The coefficient tensors obtained from microscale analysis is used for micro-macro bridging of local effect within the material. These precomputed coefficient tensors along with damage model parameters are stored as an input companion file. The external database is read in Abaqus at the beginning of the analysis, using the User-defined EXTERNAL DataBase (UEXTERNALDB) subroutine. The microscale problem is solved at each quadrature point throughout the macroscale analysis using User-defined MATerial subroutine (UMAT). The damage model parameters are evaluated within UEXTERNALDB, at the beginning of the analysis. The reduced order model problem is solved at each integration point of the macroscale numerical specimen as shown in Fig. 5. The steps employed for macroscale analysis are presented below.

5. Alleviation of spurious residual stiffness The reduced order microstructure problem constitutes a full nonlinear system used to evaluate the microstructural behavior. If low order models are used, this approach was shown to demonstrate a significant spurious post-peak residual stiffness, Crouch et al. [16] and Bogdanor et al. [18]. This spurious residual stiffness needs alleviation to ensure accurate load redistribution in the damage propagation scenario. This issue was attributed to the concept of the incompatibility of the eigenstrains as described in Furuhashi et al. [29] and Fish et al. [30]. The coefficient tensors with degraded constituent properties in the damage path represents cracking across the microstructural domain, which in turn eliminates the corresponding stiffness properties of the composite material in the appropriate loading direction. The identification of failure path and implementation of modified coefficient tensors for the ROM is described in depth by Sparks et al. [31]. For the sake of brevity, minimally required mathematical details are presented. In addition, the implementation of a generalized “m-th order” failure path approach is considered. The proposed approach needs storage of m additional sets of modified coefficient tensors, which are also computed prior to macroscale analysis. This implementation incurs insignificant additional computational cost at the macroscale. The capability of the proposed approach is demonstrated with ‘first order’ numerical verification study only. A unidirectionally reinforced composite (with the fiber oriented along Z-direction) subjected to uniaxial loading along the x-direction is considered, as illustrated in Fig. 7. The boundary conditions are chosen to be symmetric on three sides of the domain. As shown in Fig. 3, the reduced order model used in this study consist of four parts. The macroscale discretization consists

I) The numerical specimen configurations (i.e. layup, orientations, and mesh) are the input for macroscale finite element analysis using Abaqus. II) Displacement control loading is applied to the numerical specimen.

Fig. 5. Macroscale analysis with reduced order model.

7

Composite Structures 235 (2020) 111759

R. Bhattacharyya and P.K. Basu

residual stresses continued to persist due to spurious stiffness effect, ultimately culminating into failure. The progressive failure of the remaining damage parts resulted in a prolonged failure process. To remove the post-failure residual stresses, modified coefficient tensors were computed for matrix-cracking damage prior to macroscale analysis. The stress vs. strain plots in Fig. 8b shows that both irregular stress patterns and residual stresses are non-existent. This approach is extended to the model for failure of each damage part as well as in the failure of combination of several damaged parts. 6. Numerical examples In this section, calibration and validation of the proposed constitutive model are presented. For this purpose three different laminate configurations are used. Both unnotched and open-hole laminates are considered. The calibrated parameters are used in all validation simulations. 6.1. Calibration of model parameters The damage model is characterized by 11 material parameters and 4 algorithmic parameters. The algorithmic parameters are appropriately chosen to ensure performance efficiency of the model. Table 1 lists the algorithmic parameters under consideration and the nature of execution. Procedures used for evaluating these parameters are given in Appendices B and C. The constituent material parameters are calibrated for IM7/977-3 graphite epoxy composite with approximately 65% fiber volume fraction. Standardized ASTM tests in uniaxial tension [32], uniaxial compression ([0°8s ] and [90°12s ]) [33], and unnotched [+ 45°/ 45°]4s specimen were conducted at Air Force Research Laboratory (AFRL) for the purpose of calibration. Calibration of the material parameters was undertaken using experimental data from Clay and Knoth [34] and, Bogdanor and Oskay [18]. The calibrated elastic material parameters for polymer matrix and carbon fiber are listed in Table 2, 3 respectively. The calibrated constitutive model parameters for polymer matrix and carbon fiber are listed in Tables 4, 5, respectively. In all the stress vs. strain plots the simulation results are compared with 4 sets of experimental data. The unnotched specimens used for tension calibration had a length of 250 mm. On the other hand, the unnotched specimens used for compression parameter calibration had a length of 140 mm.The width of the specimen varied from 12.5 mm to 25.4 mm based on the longitudinal and transverse loading used. Thickness of each ply is 0.127 mm. In general, for all calibration cases it is observed that the simulation results show good agreement with experimental data, both in magnitude and trend. The calibration of transverse tension parameters is done using three point bend test data for [90°]8s , in place of the erroneous data for uniaxial tension shown in the Fig. 9a. In Fig. 9b, the tail end of experiment 3, is erroneous experimental artifact. The apparent overshooting of predicted response in Fig. 10 is also attributed to errors in experimental data, as conceded by the source [34]. For this reason, the compression parameters are calibrated on the basis of longitudinal compressive strength of 1680 MPa, given by manufacturer’s manual [35]. Simulation of the numerical specimens was undertaken on a parallel computing cluster comprising of 16 core 2.1 GHz Intel Xeon E5 processors with 256 GB shared memory on each compute node. Each of the simulation was performed on a single compute node using 32 cpus in a shared memory parallel configuration. [+ 45°/ 45°]4s The size of unnotched specimen is 250 mm × 38.1 mm. In the case of the numerical specimen, during shear calibration the mesh in each ply is aligned along + 45° or 45° to avoid any mesh bias in the solution. In this case as shown in Fig. 11, the stress vs. strain predictions closely follow the experimental trend, even in the softening segment.

Fig. 6. Computational implementation of the EHM framework in quasi-static failure prediction.

Fig. 7. One element model for verification.

of a single linear hexahedron. The analyses were performed using elements using both full integration and reduced integration schemes. The nodes on the positive x-face of the cube were subjected to displacement controlled loading with a maximum amplitude of 0.50 mm. In this analysis the predominant failure mode is found to be matrix cracking (i.e., the failure in matrix damage part). The stress vs. strain plot corresponding to stated loading direction is shown in Fig. 8a. This figure clearly demonstrates an irregular stress pattern beyond failure in the matrix cracking part of the reduced order model and residual stress in the macroscale element. The origin of the irregular stress pattern can be attributed to the presence of such residual stresses. Despite the failure observed in the matrix-cracking part, on further loading of the remainder of the fiber and delamination damage parts, 8

Composite Structures 235 (2020) 111759

R. Bhattacharyya and P.K. Basu

Fig. 8. Structural scale stress-strain plots under unidirectional loading. Table 1 Algorithmic parameters and their execution. Name of parameter Transition strain for uniaxial tension,

(T) tr (S) tr

Transition strain for pure shear loading, Transition strain for uniaxial compression, (C) tr

Triaxiality weighting parameter,

Table 4 Polymer matrix damage model parameters. Execution

Parameter

Value

Automated

Compression-Tension anisotropy, c Strain to initial damage, 0

0.45 0.005 0.14

Automated

Failure strain under pure shear loading,

Smooth transition in failure envelope

Shape parameter under Shape parameter under Shape parameter under Shape parameter under Shape parameter under Shape parameter under Shear stress parameter,

Value

Elastic modulus (in tension),

E (m)

(GPa)

Elastic modulus (in compression), E (m) (GPa)

Poisson’s ratio,

(m)

3.98 0.375

Value

(f ) Longitudinal elastic modulus (in tension), E33 (GPa)

252.00

Shear modulus,

(f ) G13

(GPa)

(f ) Longitudinal Poisson’s ratio, 12 (f ) Transverse Poisson’s ratio, 31

uniaxial tension loading, T uniaxial tension loading, T uniaxial compression loading, uniaxial compression loading, pure shear loading, S pure shear loading, S parm (MPa)

C

C

0.12 2.20 3.10

1.555 5.75 0.395 0.10

90.0 2.40

3.70

Parameter

(f ) Transverse elastic modulus, E11 (GPa)

(C) f

Triaxiality weighting parameter,

Table 5 Carbon fiber damage model parameters.

Table 3 Carbon fiber elastic material parameters.

(f ) Longitudinal elastic modulus (in compression), E33 (GPa)

0.38

(S) f

Failure strain under pure uniaxial compression loading,

Table 2 Polymer matrix elastic material parameters. Parameter

(T) f

Failure strain under pure uniaxial tension loading,

Automated

Parameter

Value

Compression-Tension anisotropy, c Strain to initial damage, 0

1.4481 0.0075 0.0275

Failure strain under uniaxial tension loading,

(T) f

210.50

Shape parameter under uniaxial tension loading, Shape parameter under uniaxial tension loading,

12.45

Failure strain under uniaxial compression loading,

0.5732 4.3549

T T

0.291

0.0170

(C) f

Shape parameter under uniaxial compression loading, Shape parameter under uniaxial compression loading,

13.20

C

C

0.8817 1.1686

0.206

tension, pure compression and pure shear are plotted as points on the failure envelop based on AFRL data in Table 6. Qualitatively the failure envelope appears to be very close to the failure envelope of other CFRP composites, as documented in Pinho et al. [36] and Camanho et al. [37].

Comparison of AHLS model predictions for IM7/977-3 lamina with those of AFRL are shown in Table 6. The constitutive model for polymer matrix is employed to different states of strain to verify its capability under load triaxiality. The failure envelope for the matrix material is presented in Fig. 12. The strain states at different region of the plot is schematically shown in the figure. The red cross and blue star points on the plot represent predicted ultimate strength of matrix in tension and compression regimes, respectively. The smooth transitioning is ensured by the triaxiality weighting parameter ( ) as used in Eqs. (34) and (35). It is evident from the figure that the proposed constitutive model does not need a failure criteria to predict strength. The experimental data to verify the mixed modes within failure envelope for IM7/977-3 is not available to the authors. However, the experimental data for composites under pure

6.2. Validation The unnotched specimens used for validation had a width of 25.4 mm. The length of specimen in tension was 250 mm and in compression it was 140 mm. The thickness of each ply for validation layup was 0.127 mm. In the case of open-hole specimens, the diameter of the central hole was 6.35 mm. The stress in the laminate is computed for the numerical specimens as the sum of the reactive forces at the pinned end of the specimen divided by the gross cross-sectional area of the specimen. The strain computations closely followed the method the strains 9

Composite Structures 235 (2020) 111759

R. Bhattacharyya and P.K. Basu

Fig. 9. Calibration of tension parameters.

Fig. 10. Calibration of compression parameters.

were measured in the experiments using extensometers in the case of open-hole specimens and strain-gages for unnotched specimens. For instance, two nodes in the mesh for the outer ply in each numerical

specimen were identified corresponding to the gage points of the extensometer. Then the computed strain is obtained as the change in the distance between these two nodes divided by gage length. 6.2.1. Tension In this case, all the laminates are subjected to uniaxial tensile loading along the longitudinal direction. It is apparent that the constitutive model predicts the initial stiffness quite accurately. The results for [0, 45, 90, 45]2s , [60, 0, 60]3s and [30, 60, 90, 60, 30]2s under uniaxial tension are presented in Figs. 13, 15, 17. In the case of [0, 45, 90, 45]2s unnotched and open-hole numerical specimens, the predicted stiffness is exactly same as observed in experiments. A small variation is noted in the ultimate strength prediction for open-hole specimen shown in Fig. 13a. However, the constitutive model is able to capture the nonlinearity observed near peak. Gray scale damage contours of the [0, 45, 90, 45]2s quasi-isotropic open-hole specimen for the layups up to the plane of symmetry are shown in Fig. 14. The darker shades indicate higher damage whereas the lighter ones indicate low to no damage state. The damage in the 0° plies are contributed by fiber damage, matrix cracking and delamination. On the other hand, 90° and 45° plies observed significant matrixcracking near the hole. In the case of [60, 0, 60]3s unnotched and open-hole numerical specimens, the stiffness and ultimate strength predictions are quite

Fig. 11. [+45°/ 45°]4s specimen under tension.

10

Composite Structures 235 (2020) 111759

R. Bhattacharyya and P.K. Basu

Table 6 Summary of calibrated and experimental properties of IM7/977-3 lamina. Property

Symbol

AHLS model

AFRL Mean

AFRL St.Dev

Tensile longitudinal modulus (GPa) Compressive longitudinal modulus (GPa) Tensile transverse modulus (GPa) Compressive transverse modulus (GPa) Tensile longitudinal strength (MPa) Compressive longitudinal strength (MPa) Tensile transverse strength (MPa) Compressive transverse strength (MPa) In-plane shear strength at 5% shear strain (MPa)

E1T E1C E2T E2C XT XC YT YC u F5%

164.04 137.29 8.957 8.707 2905.8 1565.6 128.0 248.56 99.8

164.0 137.0 8.98 8.69 2905.0 1569.0 44.5 248.0 99.4

4.12 0.608 0.284 0.367 146.7 44.2 2.10 11.3 2.42

1.93 1.47 3.31

1.61 1.28 3.12

0.0718 0.0586 0.124

Tensile longitudinal failure strain (%) Compressive longitudinal failure strain (%) Compressive transverse failure strain (%)

1T 1C

2C

The damage contour for the [30, 60, 90, 60, 30]2s laminate specimen shows mostly matrix-cracking in the 60° plies, as shown in Fig. 18. No significant fiber damage or delamination is observed. 6.2.2. Compression In this case, all the laminates are subjected to uniaxial compressive loading in the longitudinal direction. It is observed that the initial stiffness predicted by the proposed model closely follows the experimental data. The results for [0, 45, 90, 45]2s , [60, 0, 60]3s and [30, 60, 90, 60, 30]2s under uniaxial compression are presented in Figs. 19, 21, 23. It should be noted that in all the open-hole compression cases, the experimental data was found to be noisy, confirming the findings in [34]. The unnotched compression prediction of [0, 45, 90, 45]2s follows the experimentally observed data. However, the experimental data is beset with significant variance, as shown in 19a. The prediction for [0, 45, 90, 45]2s open-hole numerical specimen demonstrates higher strength. It is reported by ARFL (see [38]) that some of the experimental specimen failed prematurely due to buckling. Such second-order phenomenon is not considered in this study. In the [0, 45, 90, 45]2s quasi-isotropic open-hole specimen under compression the damage contour shows failure of plies in all three failure modes, as shown in Fig. 20. In the case of [60, 0, 60]3s unnotched specimen the numerical model underpredicts the strength of the laminate. The possible reasons can partly be attributed to the reduced order modeling process and partly to significant uncertainty in the experimental data, as evident from Fig. 21a. In the open-hole case, the numerical predictions shown in Fig. 21b tend to exactly follow the experimental data.

Fig. 12. Failure envelope for polymer matrix constituent.

accurate. Also, the constitutive model is able to capture the expected nonlinearity near the ultimate strength. The damage contour for the [60, 0, 60]3s open-hole laminate specimen is shown in Fig. 16. The 0° plies demonstrated fiber damage, matrix cracking and delamination. In the 60° plies fiber splitting and very small delamination are observed. The 60° ply matrix-cracking is predominant in other failure modes. In the case of [30, 60, 90, 60, 30]2s unnotched and open-hole ‘soft ply’ specimens, the stiffness predictions are also quite accurate. The numerical predictions in Fig. 17a and b closely follow the experimental data.

Fig. 13. Stress vs. strain for [0, 45, 90,

45]2s quasi-isotropic specimens under tension.

11

Composite Structures 235 (2020) 111759

R. Bhattacharyya and P.K. Basu

Fig. 14. Damage contour for [0, 45, 90,

45]2s quasi-isotropic specimen under tension.

Fig. 15. Stress vs. strain for [60, 0,

In the [60, 0, 60]3s open-hole layup fiber splitting is observed in the 60° plies. Matrix-cracking observed in 0° plies is higher than 60° plies as shown in Fig. 22. Compression stiffness of [30, 60, 90, 60, 30]2s unnotched specimen initially tended to follow the stiffness observed in the

60]3s specimens under tension.

experiment. But the reported strength for the layups were suspiciously too high for a soft layup like [30, 60, 90, 60, 30]2s . In the case of open-hole compression, the predictions agree well with experimental data except near the tail end, where the experimental data was reported to be very noisy, and the coupon specimens were loaded till 90% of the

Fig. 16. Damage contour for [60, 0,

12

60]3s specimen under tension.

Composite Structures 235 (2020) 111759

R. Bhattacharyya and P.K. Basu

Fig. 17. Stress vs. strain for [30, 60, 90,

Fig. 18. Damage contour for [30, 60, 90,

ultimate load. In the case of the [30, 60, 90, 60, 30]2s open-hole layup, damage accumulation due to matrix cracking was observed in the case of 30° and 60° angle plies, as evidenced in Fig. 24. In addition, no fiber damage or delamination is observed. This could be another reason for high strength predictions observed in Fig. 23b.

Fig. 19. Stress vs. strain for [0, 45, 90,

60,

60,

30]2s specimens under tension.

30]2s specimen under tension.

To this end, a summary of the predicted strengths is presented in Table 7 along with the simulation wall clock time on the cluster computer. The unnotched tension and compression is denoted as UNT and UNC respectively. Similarly, the open-hole specimens under tension and compression are denoted as OHT and OHC, respectively.

45]2s quasi-isotropic specimens under compression.

13

Composite Structures 235 (2020) 111759

R. Bhattacharyya and P.K. Basu

Fig. 20. Damage contour for [0, 45, 90,

45]2s quasi-isotropic specimen under compression.

Fig. 21. Stress vs. strain for [60, 0,

60]3s specimens under compression.

Fig. 22. Damage contour for [60, 0,

60]3s specimen under compression.

14

Composite Structures 235 (2020) 111759

R. Bhattacharyya and P.K. Basu

Fig. 23. Stress vs. strain for [30, 60, 90,

60,

Fig. 24. Damage contour for [30, 60, 90,

[0, 45, 90,

[60, 0,

45]2s

60]3s

[30, 60, 90,

60,

30]2s

Type

Wall clock time (hh:mm)

AHLS model (MPa)

AFRL Mean (MPa)

AFRL St. Dev (% )

UNT OHT UNC OHC

00: 03: 00: 04:

48 20 47 48

980.67 600.77 638.08 470.29

866 554 603 341

1.48 1.29 1.41 0.78

UNT OHT UNC OHC

01: 04: 01: 05:

25 04 58 47

1056.6 533.95 630.42 358.69

1005 543 765 358

1.64 1.09 2.03 0.92

UNT OHT UNC OHC

01: 01: 01: 01:

13 12 12 45

503.30 414.06 506.45 404.4

473 409 392 295

1.38 1.44 1.26 1.01

30]2s specimen under compression.

following a linear stress-strain law. The model parameters are calibrated using experimental data obtained at the laminate scale. A new triaxiality factor is defined to account for the presence of a complex strain state. Consequently, no failure criterion is needed in progressive damage modeling. The predictions of damage accumulation, stiffness, and ultimate strength of a range of laminated composite layups under tension appears to be quite satisfactory. However, modeling in compression demonstrated some discrepancy which mostly arose from the fact that the experimental data for compression parameter calibration is marred by large uncertainty. The computational efficiency is enhanced by calculating the algorithmic parameters only once at the beginning of the analysis, using a user defined subroutine UEXTERNALDB. Finally, it has been clearly demonstrated that the proposed constitutive model for polymer matrix material has the ability to capture different strain triaxiality conditions without using any phenomenological failure criteria.

Table 7 Summary of simulation wall time and predicted strengths of IM7/977-3 laminates. Layup

60,

30]2s specimens under compression.

Data availability statement Data will be made available on request.

7. Conclusions Constitutive relationships are presented for nonlinear multiscale progressive damage analysis of fiber composites. A mechanics-based formulation of constitutive relations for polymer matrix with carbon fiber reinforcement is presented. The damage evolution potential for the softening regime is defined so that the material fails consistently

CRediT authorship contribution statement Rudraprasad Bhattacharyya: Conceptualization, Investigation, Software, Writing - original draft. Prodyot K. Basu: Supervision, Writing - review & editing.

15

Composite Structures 235 (2020) 111759

R. Bhattacharyya and P.K. Basu

Declaration of Competing Interest

Directorate of the Air Force Research Laboratory (Contract No: GS04T09DBC0017 through Engility Corporation and FA8650-12-D3212/0016 through UDRI). Valuable comments by Dr. Stephen B. Clay of AFRL towards the development of the constitutive model is acknowledged. In addition, Dr. Caglar Oskay’s intellectual support and access to HPC resources to the first author is acknowledged.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgement Partial financial support was provided by the Aerospace Systems Appendix A. Derivation of the triaxiality factor

The concept of triaxiality factor is to control the constitutive model parameters to capture the nonlinear shear behavior to brittle uniaxial behavior. The old definition of triaxiality factor was given by Eq. (30). In this work, the range of triaxiality is varied from purely uniaxial to pure shear loading. However, the numerical range of triaxiality factor remains the same. Keeping in view the two strain-states of interest, the expressions of new triaxiality factor is derived. A.1. Pure uniaxial loading Assuming an isotropic material, under pure uniaxial tension loading, the principal strain vector is written as: { shear strain max and maximum absolute principal strain max are computed as Eqs. (A.1) and (A.2) max

2 max

( 2

=

)

=

(1 + ) 2

} . The maximum

(A.1) (A.2)

=

Substituting the expressions in Eq. (30) the triaxiality factor is computed as:

kb =

(1 + ) (1 + ) 2

=

+

2(1 + ) (3 + )

(A.3)

In ensuring kb = 0 under pure uniaxial loading, the expression takes the form of Eq. (A.4).

kb =

max

(

max

2

+

max

2(1 + ) (3 + )

)

(A.4)

This treatment is verified to be true both for uniaxial tensile and compressive loading. A.2. Pure shear loading In order to ensure kb = 1 under pure shear loading, Eq. (A.4) needs to be factorized by x. Under pure shear loading the maximum shear strain and maximum absolute principal strain max are computed as in Eqs. (A.5) and (A.6) max

2 max

( 2

=

)

=

max

(A.5) (A.6)

=

Substituting the maximum strains in Eq. (A.4) with a factor x gives the following expression:

1 =x

2 +

=x

3+

2(1 + ) (3 + ) 2 3+

2

Hence, the factor x is:

3+ 1

x=

Consequently, The final expression for triaxiality factor takes the form of Eq. (A.7).

kb =

(3 + ) (1 )

(

max max

2

+

max

)

2(1 + ) (3 + )

(A.7)

Appendix B. Evaluation of damage equivalent strain parameters for matrix The damage equivalent strain at each load increment of loading is obtained using Eq. (25). In order to calculate the damage equivalent strain parameters (3 for pure uniaxial and 3 for pure shear), the corresponding two strain state vectors are needed. First, the expressions corresponding to pure uniaxial tension and pure shear strain states are obtained. The strain vector for pure uniaxial tension 16

Composite Structures 235 (2020) 111759

R. Bhattacharyya and P.K. Basu

and pure shear loading are given by Eqs. (B.1) and (B.2), respectively. (N) (·) (N) (·)

=

0

0 0

(S) (·)

0 (N) (·)

(N) (·)

0

0

(S) (·)

0

(S) (·)

0

0

0

0

0

=

0 (B.1)

(B.2)

The principal strain vectors for pure uniaxial tension and pure shear loading are obtained as Eqs. (B.3) and (B.4), respectively. (N)

=

{

(S)

=

{

(N) (·)

(N) (·)

(S) (·)

(S) (·)

0

(N) (·)

}

(B.3)

}

(B.4)

B.1. Pure uniaxial loading Substituting the expression for principal strain vector in Eq. (B.3) into the strain tensor of Eq. (25) the expression of damage equivalent strain parameter can be written as Eq. (B.5)

(N) (·)

1 h1 2

=

(N) (·)

The final expression of (N) (·)

=

(N) (·)

(N) (·)

h2

(N) (·)

+ 2µ + 2µ

h2

(N) (·)

h3

(N) (·)

(B.5)

after algebraic operations becomes Eq. (B.6).

(N) (·)

E 2(1 + )(1

h3

(N) (·)

h1

+ 2µ

)(h12 +

(1

2 )

2h 2 2

+

2h 2 ) 3

2 2 (h1 h2 + h1 h3

h2 h3)

(B.6)

Following the logic of Eq. (27), h1 = 1 and h2 = h3 = c . Consequently, the Eq. (B.6) finally takes the form of Eq. (36). B.2. Pure shear loading The development of damage equivalent strain parameters for shear loading follows the same methodology. In this case the damage equivalent strain parameter is expressed as Eq. (B.7). (S) (·)

=

h1

(S) (·)

0 h3

(S) (·)

The final expression of S (·)

=

S (·)

h1

+ 2µ

1 2

E 2(1 + )(1

+ 2µ

0 + 2µ

h3

(S) (·)

(B.7)

after algebraic operations becomes Eq. (B.8).

(S) (·)

)(h12 + h32 )

(1

2 )

(S) (·)

2 h1 h3

(B.8)

Following the logic of Eq. (27), h1 = c and h3 = 1. Hence, the Eq. (B.8) finally takes the form of Eq. (37). The damage equivalent strain parameters at damage threshold or initiation and failure are computed using experimental data. Evaluation of transition damage equivalent strain for pure normal and pure shear loading can be computed using the analytical expressions. The transition damage equivalent strain is physically conceived as a strain at which (or at vicinity of which) the stress-strain diagram reaches a peak at the onset of softening. Mathematically, this point can be obtained by making the gradient of stress state equals to zero (Eq. (B.9)).

d d

= 1

C:

d d

d C: d

=0

(B.9)

Evaluation of for pure normal loading is straightforward as there is a well-defined peak in the stress vs. strain plot. This parameter ( obtained by solving Eq. (B.10): (N) tr

N(m)

=E

1

(

(N) N tr

(N) tr ) 2

+ atan(

N)

[1 + (

N

(N) tr

0

N)

2]

(N) tr )

is

=0 (B.10)

However, the above methodology cannot be directly followed for pure shear loading as there is no such peak in the hardening regime due to presence of high ductility. The transition parameter, tr(S) for this case is obtained from shear stress parameter ( parm ), which is a parameter of this model. Eq. (B.11) is iteratively solved to obtain the transition parameter as:

17

Composite Structures 235 (2020) 111759

R. Bhattacharyya and P.K. Basu

Fig. B.25. Evaluation of transition damage equivalent strain parameters.

S(m )

=

atan(

S

(S) tr

S)

0

+ atan( S)

+ atan( S)

2

parm (1

+

E

+ )

S

1=0

(S) tr

(B.11)

Evaluation of the transition damage equivalent strain parameters using Eqs. (B.10) and (B.11) are presented in Fig. B.25a and b, respectively. Appendix C. Evaluation of damage equivalent strain parameters for fiber The model parameters that describe damage evolution and transitioning from the damage induced strain hardening regime to softening is obtained based on the analysis of the evolution equations under normal loading along the direction of the fiber length (assumed as the 3-direction in what follows). The strain tensor at initiation, transition or at final failure stages under uniaxial loading is written as: 31 (·)

0

(·)

0 0

=

31

0 0

(·)

0

(C.1)

(·)

The principal strains for the pure uniaxial loading are therefore:

={

31

(·)

31

(C.2)

(·) }

(·)

Substituting Eqs. (25)–(27) into the expression for principal strain vector in Eq. (C.2), the expression of damage equivalent strain parameter can be written as:

(·)

=

1 [ 2

31 c

31 c

(·)

(·)

(·) ]

31 c

C11 C12 C13 C12 C11 C13 C13 C13 C33

31 c

(·) (·)

(C.3)

(·)

Eq. (C.3) can further be simplified by matrix operations, as below.

(·)

= =

(·)

(·)

1 [ 2 1 2

31 c

31 c

1]

31 cC11

31 cC12

31 cC12

31 cC11

+ C13 + C13 31 cC12 + C33

31 cC13

2

2 2 31 c C11

+2

2 2 31 c C12

4

31 cC13

After rearranging, the final expression of (·)

=

(·)

C33 +

2 2c 2 31 (C11

+ C12) 2

4c

(·)

31 C13

+ C33

(C.4)

is expressed in the following form as in Eq. (42):

=

(·)

(f ) N

The computation of transition damage equivalent strain for carbon fiber can be obtained from Eq. (B.9), similar to polymer matrix. Evaluation of for normal loading is straightforward as there is a well-defined peak in the stress vs. strain plot. The gradient is computed from the ascending part of the curve (i.e., arctangent law). Considering the fiber subjected to uniaxial loading along along its length, the final expression for Eq. (B.9) in terms of tr becomes Eq. (C.5): tr

(f )

= [C33

2

31 C13]

1

(

tr )

tr

[ 2 + atan( )][1 + (

tr

0

) 2]

18

=0

(C.5)

Composite Structures 235 (2020) 111759

R. Bhattacharyya and P.K. Basu

Appendix D. Supplementary data Supplementary data associated with this article can be found, in the online version, at https://doi.org/10.1016/j.compstruct.2019.111759.

loading. J Compos Mater 2017;51(10):1455–72. [19] Soutis C, Turkmen D. Moisture and temperature effects of the compressive failure of cfrp unidirectional laminates. J Compos Mater 1997;31(8):832–49. https://doi.org/ 10.1177/002199839703100805. [20] Gilat A, Goldberg RK, Roberts GD. Strain rate sensitivity of epoxy resin in tensile and shear loading. J Aerospace Eng. [21] Littell JD, Ruggeri CR, Goldberg RK, Roberts GD, Arnold WA, Binienda WK. Measurement of epoxy resin tension, compression, and shear stress–strain curves over a wide range of strain rates using small test specimens. J Aerospace Eng 2008:162–73. [22] Park IK, Park KJ, Kim SJ. Rate-dependent damage model for polymeric composites under in-plane shear dynamic loading. Comput Mater Sci 2015;96:506–19. [23] Salavatian M, Smith LV. An investigation of matrix damage in composite laminates using continuum damage mechanics. Compos Struct 2015;131:565–73. [24] Van Paepegem W, De Baere I, Degrieck JJ. Modelling the nonlinear shear stressstrain response of glass fibre-reinforced composites. Part I: experimental results. Compos Sci Technol 2006;66(10):1455–64. https://doi.org/10.1016/j. compscitech.2005.04.014. URL:http://www.sciencedirect.com/science/article/pii/ S0266353805000965. [25] Lomakin EV, Fedulov BN. Nonlinear anisotropic elasticity for laminate composites. Meccanica 2015;50(6):1527–35. https://doi.org/10.1007/s11012-015-0104-5. [26] Van Paepegem W, De Baere I, Degrieck JJ. Modelling the nonlinear shear stressstrain response of glass fibre-reinforced composites. Part II: model development and finite element simulations. Compos Sci Technol 2006;66(10):1465–78. https://doi. org/10.1016/j.compscitech.2005.04.018. URL:http://www.sciencedirect.com/ science/article/pii/S026635380500117X. [27] Fedulov B, Fedorenko A, Safonov A, Lomakin E. Nonlinear shear behavior and failure of composite materials under plane stress conditions. Acta Mech 2017;228(6):2033–40. https://doi.org/10.1007/s00707-017-1817-4. [28] Feyel F, Chaboche JL. FE2 multiscale approach for modelling the elastoviscoplastic behavior of long fiber SiC/Ti composite materials. Comput Methods Appl Mech Eng 2000;183:309–30. [29] Furuhashi R, Mura T. On the equivalent inclusion method and impotent eigenstrains. J Elasticity 1979;9(3):263–70. [30] Fish J, Filonova V, Yuan Z. Hybrid impotent-incompatible eigenstrain base homogenization. Int J Numer Methods Eng 2013;95(1):1–32. [31] Sparks PA, Oskay C. The method of failure paths for reduced-order computational homogenization. Int J Multiscale Comput Eng 2016;14:515–34. https://doi.org/10. 1615/IntJMultCompEng. 2016018702. [32] ASTM Standard 3039. Standard test method for tensile properties of polymer matrix composite materials. West Conshohocken, PA: ASTM International; 2008. [33] ASTM Standard D3410. Standard test method for compressive properties of polymer matrix composite materials with unsupported gage section by shear loading. West Conshohocken, PA: ASTM International; 2008. [34] Clay SB, Knoth PM. Experimental results of quasi-static testing for calibration and validation of composite progressive damage analysis methods. J Compos Mater 2017;51(10):1333–53. [35] Cytec Engineered Materials. CYCOM 977-3 epoxy resin system technical data sheet; 2012. URL:https://www.cytec.com/sites/default/files/datasheets/CYCOM_977_3. pdf. [36] Pinho ST, Davila CG, Camanho PP, Iannucci L, Robinson P. Failure models and criteria for FRP under in-plane or three-dimensional stress states including shear non-linearity, techreport WU 23-064-30-22. NASA; 2005. URL:https://ntrs.nasa. gov/archive/nasa/casi.ntrs.nasa.gov/20050110223.pdf. [37] Camanho PP, Arteiro A, Melro AR, Catalanotti G, Vogler M. Three-dimensional invariant-based failure criteria for fibre-reinforced composites. Int J Solids Struct 2015;55:92–107. https://doi.org/10.1016/j.ijsolstr.2014.03.038 Special Issue Computational and Experimental Mechanics of Advanced Materials A workshop held at King Abdullah University of Science and Technology Jeddah, Kingdom of Saudi Arabia July 1–3, 2013. URL:http://www.sciencedirect.com/science/article/ pii/S0020768314001449. [38] Engelstad SP, Clay SB. Comparison of composite damage growth tools for static behavior of notched composite laminates. J Compos Mater 2017;51(10):1493–524. https://doi.org/10.1177/0021998316675945.

References [1] Department of Defense, Composite Materials Handbook, F Edition, vol. 2. Polymer Matrix Composites Materials Properties. SAE International; 2002. [2] NASA. Advanced composites partnership [Tech. rep.]. NASA; 2014. [3] Liu X, Furrer D, Kosters J, Holmes J. Vision 2040: A roadmap for integrated, multiscale modeling and simulation of materials and systems [Technical Report NASA/ CR-2018-219771, E-19477, GRC-E-DAA-TN52454]. NASA; 2018. [4] Pinho ST, Iannucci L, Robinson P. Physically-based failure models and criteria for laminated fibre-reinforced composites with emphasis on fibre kinking: part I: development. Compos Part A 2006;37(1):63–73. https://doi.org/10.1016/j. compositesa.2005.04.016. URL:http://www.sciencedirect.com/science/article/pii/ S1359835X05002198. [5] Pinho ST, Iannucci L, Robinson P. Physically based failure models and criteria for laminated fibre-reinforced composites with emphasis on fibre kinking. Part II: FE implementation. Compos Part A 2006;37(5):766–77. https://doi.org/10.1016/j. compositesa.2005.06.008. [6] Ridha M, Wang CH, Chen BY, Tay TE. Modelling complex progressive failure in notched composite laminates with varying sizes and stacking sequences. Compos Part A 2014;58:16–23. https://doi.org/10.1016/j.compositesa.2013.11.012. URL:http://www.sciencedirect.com/science/article/pii/S1359835X13003217. [7] Leone FA, Dávila CG, Girolamo D. Progressive damage analysis as a design tool for composite bonded joints. Compos Part B 2015;77:474–83. https://doi.org/10.1016/ j.compositesb.2015.03.046. URL:http://www.sciencedirect.com/science/article/ pii/S1359836815001651. [8] Su ZC, Tay TE, Ridha M, Chen BY. Progressive damage modeling of open-hole composite laminates under compression. Compos Struct 2015;122:507–17. https:// doi.org/10.1016/j.compstruct.2014.12.022. URL:http://www.sciencedirect.com/ science/article/pii/S0263822314006783. [9] Higuchi R, Okabe T, Nagashima T. Numerical simulation of progressive damage and failure in composite laminates using XFEM/CZM coupled approach. Compos Part A 2017;95:197–207. https://doi.org/10.1016/j.compositesa.2016.12.026. URL:http://www.sciencedirect.com/science/article/pii/S1359835X16304547. [10] van Dongen B, van Oostrum A, Zarouchas D. A blended continuum damage and fracture mechanics method for progressive damage analysis of composite structures using XFEM. Compos Struct 2018;184:512–22. https://doi.org/10.1016/j. compstruct.2017.10.007. URL:http://www.sciencedirect.com/science/article/pii/ S0263822317318585. [11] FLaurin F, Carrère N, Maire J-F. A multiscale progressive failure approach for composite laminates based on thermodynamical viscoelastic and damage models. Compos Part A 2007;38(1):198–209. https://doi.org/10.1016/j.compositesa.2006. 01.018. URL:http://www.sciencedirect.com/science/article/pii/ S1359835X0600025X. [12] Pineda EJ, Bednarcyk BA, Arnold SM. Validated progressive damage analysis of simple polymer matrix composite laminates exhibiting matrix microdamage: comparing macromechanics and micromechanics. Compos Sci Technol 2016;133:184–91. https://doi.org/10.1016/j.compscitech.2016.07.018. URL:http://www.sciencedirect.com/science/article/pii/S0266353816307692. [13] Massarwa E, Aboudi J, Haj-Ali R. A multiscale progressive damage analysis for laminated composite structures using the parametric HFGMC micromechanics. Compos Struct 2018;188:159–72. https://doi.org/10.1016/j.compstruct.2017.11. 089. URL:http://www.sciencedirect.com/science/article/pii/ S026382231730702X. [14] Fish J. Practical multiscaling. John Wiley & Sons; 2013. [15] Oskay C, Fish J. Eigendeformation-based reduced order homogenization for failure analysis of heterogeneous materials. Comput Methods Appl Mech Eng 2007;196:1216–43. [16] Crouch RD, Oskay C. Symmetric mesomechanical model for failure analysis of heterogeneous materials. Int J Multiscale Comput Eng 2010;8:447–61. [17] Crouch RD, Clay SB, Oskay C. Experimental and computational investigation of progressive damage accumulation in CFRP composites. Compos Part B 2013;48:59–67. [18] Bogdanor MJ, Oskay C. Prediction of progressive damage and strength of IM7/9773 composites using the Eigendeformation-based homogenization approach: static

19