On Micro-Buckling of Unidirectional Fiber-Reinforced Composites by Means of Computational Micromechanics

Abstract

Micro-buckling of unidirectional fiber-reinforced composites is investigated in this paper by means of an explicit representation of a geometrically imperfect fiber within the context of kinematical and material non-linear behavior. Two types of fiber imperfections are considered: a helicoidal shape, identified as 3D imperfection; and a sinusoidal plane shape (2D imperfection). Both imperfection models are characterized by a maximum misalignment angle of the fiber with respect to the ideal or perfect configuration, as is usually considered in this field. A total of 816 cases were computed in terms of imperfection type (either 2D or 3D), fiber volume fraction, fiber arrangement (square or hexagonal array), orientation for 2D models, matrix yield stress, and misalignment angle. Two load cases, with constrained and unconstrained transverse strain, were considered. Assuming periodic boundary conditions, homogenization was carried out to obtain macroscopic stresses. Numerical results are compared with an analytical model available in the literature. The results show a high imperfection-sensitivity for small misalignment angles; on the other hand, the type of imperfection and the fiber arrangement do not have a large influence on the results. In addition, it was found that this problem is governed by fiber volume fraction and matrix yield stress only for small imperfections, whereas for large misalignment angles, a change in fiber volume fraction produces small changes in micro-buckling stress.

Keywords:
Composites; micromechanics; fiber misalignment; micro-buckling

1 INTRODUCTION

In the context of fiber-reinforced composite materials, the term fiber micro-buckling refers to the buckling of fibers involving transverse displacements under compression in the direction of the fiber. By analogy with a structural behavior, micro-buckling is frequently modeled as the buckling of a column which is laterally supported on an elastic matrix. This problem has attracted a number of researchers, starting from the pioneer work of Rosen (1965). There are two general ways in which such instability can occur: either due to elastic buckling of the fiber involving deformations of the matrix (a problem which is usually called micro-buckling), or by plastic deformations (called kinking), as reported by Budiansky and Fleck (1993). Comparison between both modes has been made by Sun and Tsai (2001) among others: both are seen to be different problems and involve different assumptions; however, these authors concluded that the loads (stress fields) found at buckling or kinking are not very different.

Micro-buckling may occur in mainly two modes: in a periodic (or "in-phase") mode there are periodic deformations of the fibers and shear of the matrix but no significant transverse deformation occurs. Micro-buckling can also occur in a non-periodic (or "out-of-phase") mode, in which shear is negligible and buckling is dominated by transverse deformations, much in the form of a beam on an elastic foundation. Both cases were discussed by Rosen (1965) and were reviewed in the classical text by Jones (1975). Regarding stability of the post-critical path, Maewal (1981) anticipated a stable post-buckling behavior for the in-phase problem. On the other hand, Yurgartis and Sternstein (1994) performed tests that showed the detrimental influence of fiber misalignment on the micro-buckling process, a behavior that is typical of unstable postbuckling.

Tomblin et al. (1997) investigated imperfection-sensitivity of fiber micro-buckling as a consequence of fiber misalignment under compression in the direction of the fiber within the framework of the general theory of elastic stability. Material nonlinearity was introduced in the shear constitutive equations by means of a hyperbolic relation, and the adopted RVE focused on in-phase modes. This led to the identification of an unstable symmetric bifurcation behavior for a RVE, in which a 2/3 power law was found to characterize the imperfection-sensitivity caused by fiber misalignment. Such deterministic theoretical results indicated the nature of the expected unstable behavior and were followed by probabilistic studies by Tomblin and Barbero (1997) and Barbero (1998).

The importance of stacking sequence on micro-buckling was highlighted by Drapier et al. (1996)with reference to bending tests which were performed to identify limit states under compression. Rather than considering a UC, these authors took into account a two-dimensional model of a homogenized laminated composite. For the same periodic configuration of fiber misalignment (under an in-phase mode), a nonlinear Finite Element model of a unidirectional carbon-fiber and epoxy-matrix composite was used (Drapier et al., 1998). The results under compression with misalignment showed significant loss of stress carrying capacity, exhibiting high imperfection-sensitivity. This was the first numerical study to quantify the effect of geometric imperfections on micro-buckling.

Consideration of waviness in other composites, such as biaxial and triaxial textile composites, was recently addressed by Mallikarachchi et al. (2013), Kueh (2013, 2014), and Rasin et al. (2016).

All models reported in the literature were based on a two-dimensional idealization of fiber waviness, where as in reality the fiber may adopt a 3D helicoidal shape rather than a 2D sinusoidal one. Nevertheless, the 2D idealization is often used to calculate imperfection sensitivity of fiber micro-buckling because it is mathematically tractable. Therefore, a pressing question is whether or not a 2D approximation is sufficiently accurate with respect to a 3D one, and thus the motivation for this workin which a more detailed model and refined results are discussed.

2 MODEL

2.1 Framework of Analysis

At a unit cell level, a composite is modeled in this work by means of computational micro-mechanics (see, for example, Zohdi and Wriggers, 2008);in which there are macro and micro scales. This methodology has been employed to model a wide range of heterogeneous materials, such as particle-reinforced composites (Li and Wongsto, 2004); and fiber-reinforced composites (Car et al., 2002). Matrix and fiber are modeled as two separate materials at the micro scale; whereas a single homogenous material (which is assumed to have an equivalent behavior to the heterogeneous fiber-matrix material) is taken into account at the macro scale. In a periodic heterogeneous material, this UC is employed to reconstruct the composite material by means of a repetitive pattern, as shown in Figure 1.

Figure 1:
Unit cell Example of a two-dimensional periodic material.

In computational micro-mechanics, a UC is considered with given boundary conditions and under a specific load configuration. With the solution of the stress field at the UC level, usually known as microscopic stresses, a post-process follows to evaluate the stress field in an equivalent homogeneous material, usually known as macroscopic stresses. This process at the micro level is carried out in this work using a Finite Element approximation. The model employed in this work is discussed in this section, including UC geometry, boundary conditions, and constitutive materials. The stability analysis is performed by means of the general purpose Finite Element code ABAQUS (2009).

2.2 Unit Cell Geometry

A UC shown in Figure 2a was used in this work to represent the microstructure in the periodic composite together with the geometric deviations with respect to an ideal or perfect configuration. The UC is constructed by means of the scan of a transverse section (see Figure 2b) by introducing displacements without rotation along the curve that defines the fiber along the imperfect location (imperfection lines, shown in Figure 3).Finite Element models were investigated by assuming hexagonal (designated as Hx) and square (Sq) fiber arrangements, which are shown in Figure 2b. This allows representation of a three-dimensional periodic microstructure, leading to a geometry that can be easily meshed.

Figure 2:
(a): Unit cell. (b): Transverse sections (fiber arrangements).

Figure 3:
Definition of imperfection. (a): 2D model; (b) 3D model.

For a given fiber radius Rf, the dimensions of the UC were calculated based on fiber volume fraction V_f and the ratio between fiber length Lf and diameter Df, i.e. β= Lf/Df. Since a first order computational micromechanics technique is employed in this work, absolute UC dimension values, such as fiber radius Rf, would not affect results. Higher order homogenization theories would be necessary to take into account such absolute UC dimensions, as explained by van Dijk (2016) among others. The angle θ is equal to 60º in Hx, and 90º in the Sq configuration (see Figure 2). The value of b in Figure 2 is given by

(1)

2.3 Fiber Imperfections

Two types of imperfections with respect to a straight fiber were considered in this work: (i) Three-dimensional (3D) deviations, with the fiber having a helicoidal shape; (ii) Two-dimensional (2D) deviations following a sinusoidal shape. In their initial positions, the fibers are assumed to be in phase. As shown in Figure 3, the 3D imperfection is contained in a cylinder whereas the 2D imperfection lies in a plane.

With reference to Figure 3, perfect alignment would be given by a fiber placed in the direction of axis x₁ ; a 2D imperfection has been illustrated in the plane containing axis x₁ and has a given angle δ measured with respect to the plane x₁ -x₂ ; in other words, δ defines the orientation of the plane where the 2D imperfection develops. Following the usual definitions in the micro-buckling field, both imperfections have a maximum angle of misalignment (α) of the fiber with respect to axis x₁ .

The difference between angles δ and α should be emphasized. The orientation of the plane which contains the imperfect fiber in the 2D model is δ; whereas α is the maximum angle of the fiber with respect to axis x₁ . Thus, angle δ has relevance in a 2D model, but α is a relevant parameter in both 2D and 3D models.

The Cartesian equations of the 2D imperfection line are given by

(2)

The equations for a 3D imperfection line are

(3)

where r_c is the radius of the cylinder which contains the helicoid. This value may be obtained from the expression

(4)

2.4 Periodic Boundary Conditions

Periodic boundary conditions (PBC) were used to represent a periodic microstructure. PBC are described at large in the literature on computational micro-mechanics, such as Li and Wongsto (2004), Sharma et al. (2014),among others; they were also used by Kueh and Pellegrino (2007) and by Mallikarachchi et al. (2011) by means of 1D beam elements for 2D structures.

To model the micro-structure in a periodic material it is possible to employ the concept of periodicity vectors (Oller et al., 2005; Car et al., 2002). Following the nomenclature adopted in Zahr-Viñuela and Pérez-Castellanos (2011), two points in a microstructure are identified as "corresponding points" if the position of one of them may be obtained as the position of the other one plus a linear combination of the periodicity vectors using integer coefficients. To illustrate the concept, periodicity vectors P₁ and P₂ are shown in Figure 4. The points in pairs:(C₀; C₁), (C₀; C₂) and (C₀; C₃) are corresponding points.

Figure 4:
Corresponding points.

Three vectors of periodicity, shown in Figure 5, were used in this work.

(5)

where i, j, k are unit vectors in the coordinate directions x₁ , x_{2 , and} x_{3 .}

Figure 5:
Periodicity vectors used in this work.

The boundary conditions are relations involving the forces and displacements at the boundary of the cell (Li and Wongsto, 2004). If the traction at a boundary point and its corresponding point are t⁺ and t^- respectively, then the following condition should be satisfied at all boundary pairs of points

(6)

Assuming that the displacements at a boundary point are written as u⁺ and at its corresponding point as u^-, then the condition

(7)

applies at all boundary points, where ε⁰ is the macroscopic strain tensor; and P is a periodicity vector (or a linear combination of them) which satisfies the condition

(8)

where X⁺ and X^- are the coordinates of each node of the considered couple, for which the displacements are u⁺ and u^-.

For a UC modeled by Finite Elements, it is only necessary to specify the conditions (7) for the displacements at the boundary. The conditions (6) for the boundary tractions are automatically satisfied because a displacement-based variational Finite Element formulation is employed, as explained by Li and Wongsto (2004).

There are several ways to implement the conditions (7) in practice, including Lagrange multipliers or a penalty formulation. In this work, periodicity was implemented by means of multipoint linear constraints (using the *EQUATION command in ABAQUS). Basically, following(7), the scalar equation

(9)

(for i=1, 2, 3) holds, where U_i ^j are the displacement components in direction i of the additional node j that has been included as a control node. Three control nodes have been selected in this work, and a boundary condition is assigned to them in terms of the components of the macroscopic strain tensor which is to be imposed on the UC:

(10)

To avoid problems with units, equations (9) and (10) are assumed to include a unit coefficient to homogenize units. Equations (7) are thus implemented by use of equation (9) and boundary conditions (10).

Some further details need to be considered. When implemented using multipoint linear constraints in a program like ABAQUS, the first degree-of-freedom involved in equation (9) is eliminated as an unknown from the problem and cannot be further employed in new restrictions or in other boundary conditions. For this reason, the equations that are included in the programming should be selected to achieve an effective enforcement of all necessary relationship at the UC boundary. Following this procedure it is possible to select the conditions to be employed and avoid problems with missing degrees-of-freedom.

Boundary faces, edge lines, and vertices, are identified in Figure 6.

Figure 6:
Identification of (a) faces, (b) edges, and (c) vertices.

All equations to be implemented have the form of equation (9), but for each pair of nodes itis necessary to specify the associated vector P. As an example, consider a couple of faces, such as R and L in Figure 6, and corresponding nodes on each face. Equation (9) establishes the condition for this pair of nodes, and for any other pair of corresponding nodes on faces R and L with vector of periodicity P1:

(11)

A summary of the P vectors to be used with each pair of nodes on faces, edges, and vertices is given in Tables 1 to 3 for the first and second degree of freedom (DOF) in equation (11). This procedure allows implementation of any macroscopic strain.

Table 1, 2 and 3
: P vectors of different pairs of faces, edges and vertices.

1st DOF	2nd DOF	P		1st DOF	2nd DOF	P		1st DOF	2nd DOF	P
R	L	P1		V2	V1	P1		E3	E1	P2
S	I	P2		V5	V1	P2		E11	E1	P3
F	P	P3		V3	V1	P3		E12	E1	P2+P3
				V4	V1	P1 + P3		E8	E5	P1
				V7	V1	P2 + P3		E7	E5	P2
				V6	V1	P1 + P2		E10	E5	P1 +P2
				V8	V1	P1 +P2 + P3		E4	E2	P1
								E6	E2	P3
								E9	E2	P1+P3

2.5 Load Cases

Under compression in direction x₁ , a single lamina can freely deform in the transverse directions. However, if it is used as part of a laminate, then it cannot deform transversely in the same way because other laminae produce an effect of transverse constraint. The question arises regarding how this may influence the micro-buckling of a fiber which is part of a laminate? To distinguish between possible situations, two loading cases are considered in this work: one with constrained transverse strain and another one without such constraint. Notice that the first case is not fully representative of the strain state of a lamina in a laminate, but it may be seen as a limiting condition.

Both loading cases are deformations without macroscopic distortions, i.e. the macroscopic strain tensor has zero off-diagonal components and a compressive component 11. In one case (which is identified as case A) the transverse strain components are not specified, so that the UC can expand in the transverse directions while being under 11 compression. In the other load case (identified as case B) the transverse strain components are zero, with the consequence that no lateral expansion occurs under axial compression.

Implementation of load cases A and B described above is done through boundary conditions, eq. (10), applied on displacements at control nodes since they are directly related to macro strain components. In load case B, all displacements are specified; this means that all macroscopic strain components ε⁰ _ij are imposed. For load case A, boundary conditions related to ε⁰ ₁₁ diagonal component and zero off-diagonal components are imposed; whereas the displacements of control nodes related to transverse strains ε⁰ ₂₂ and ε⁰ ₃₃ are kept as degree of freedoms in the model, and they can be obtained from the numerical solution of the UC.

2.6 Fiber and Matrix Materials

In order to investigate imperfection sensitivity, this model includes two sources of non-linearity: The epoxy matrix is modeled as a non-linear material, whereas the strain-displacement relations are geometrically non-linear. The nonlinear equilibrium paths are followed using the Riks algorithm in ABAQUS (2009).

To illustrate the micro-buckling behavior, a composite made of glass fibers and epoxy matrix was considered having R_f = 3.5 x10^-6 m. Both constituent materials were assumed to be isotropic: the fiber was modeled as linear elastic with modulus E_f = 84 GPa, and Poisson's ratio ν_f = 0.22, whereas the matrix was modeled as an elasto-perfectly plastic material with E_m= 4 GPa and ν_m = 0.38. Two yield stresses for the matrix material were considered with values given by σ_y = 48.26 MPa and σ_y =100 MPa. An associative flow rule was used together with von Mises yield surface.

2.7 Post-Processing and Finite Element Mesh

The macroscopic Cauchy stress tensor, σ^k _ij , is computed, as in many other works, as

(12)

where σ^k _ij is the ij component of the microscopic Cauchy stress tensor at Gauss point k in the Finite Element mesh that covers the UC; V_k is the integration weight (in terms of volume associated with the Gauss point k) for a mesh with N Gauss points and, V is the current volume of the UC.

Approximately 30,000 elements (C3D8 in the ABAQUS library, a solid 8-node linear brick) were used in the Finite Element mesh. This element does not have volumetric locking when used in plasticity problems. Meshes with 200,000 elements were also used in some configurations to test convergence. The maximum stress changed less than 0.7% respect to values reported in the following Section.

3 RESULTS AND DISCUSSION

Under load condition A, 384 cases were solved for this work, where as 432 cases were computed for load condition B. The cases cover changes in the variables of interest, which include the type of imperfection (either 2D or 3D); orientation angle δ of the 2D imperfection; configuration of fiber arrangement; matrix yield stress σ_y; fiber volume fraction V_f , and misalignment angle α. Cases with V_f = 10%, 30%, 50% and 70% were investigated, together with values of angle α between 0.01º and 20º and β = Lf/Df between 12.5 and 200 to have a wide perspective of the phenomenon covering results reported in experimental tests; as a reference value, Jochum and Grandidier (2004) measured values of β between 20 and 50 with misalignment angles up to 10º.

To facilitate the presentation of results, each case was identified with a group number ranging from 1 to 24, a letter to identify the load case, the misalignment angle, and the yield stress of the matrix. Group numbers are given for a V_f value, a type of imperfection, a specific fiber arrangement, and an angle δ (in cases of 2D imperfection). The codification is shown on Tables 4 to 7. For example, Group 4 includes cases with V_f = 70%, Sq arrangement, 2D imperfection, and δ= 0.

Table 4:
Limit stresses in [GPa] for unconstrained load case A and σ_y = 48.26 MPa.

Nº	Vf [%]	Imperfection	Array	delta	Misalignement angle a
					0.01º	0.1º	1º	2.5º	5º	10º	15º	20º
1	70	3D	Hx	-	3.706	2.591	1.023	0.558	0.318	0.177	0.130	0.106
2			Sq	-	3.960	2.674	1.011	0.552	0.321	-	-	-
3		2D	Hx	0º	3.752	2.639	1.158	0.592	0.372	-	-	-
4			Sq	0º	3.987	2.677	1.017	0.561	0.335	-	-	-
5			Hx	30º	3.624	2.636	1.087	0.575	0.336	-	-	-
6			Sq	45º	3.965	2.794	1.179	0.627	0.358	-	-	-
7	50	3D	Hx	-	2.303	1.702	0.813	0.479	0.286	0.162	0.116	0.093
8			Sq	-	2.329	1.681	0.809	0.479	0.287	0.162	0.117	0.094
9		2D	Hx	0º	2.346	1.744	0.867	0.519	0.318	0.207	-	-
10			Sq	0º	2.337	1.724	0.756	0.487	0.304	0.195	-	-
11			Hx	30º	2.333	1.722	0.824	0.495	0.306	0.201	-	-
12			Sq	45º	2.379	1.787	0.878	0.516	0.312	0.241	-	-
13	30	3D	Hx	-	1.480	1.146	0.626	0.398	0.256	0.152	0.110	0.088
14			Sq	-	1.505	1.157	0.624	0.397	0.255	0.151	0.110	0.088
15		2D	Hx	0º	1.525	1.175	0.656	0.426	0.284	0.184	0.149	0.128
16			Sq	0º	1.515	1.150	0.643	0.417	0.280	0.179	0.142	0.121
17			Hx	30º	1.521	1.161	0.651	0.419	0.281	0.181	0.147	0.128
18			Sq	45º	1.529	1.177	0.660	0.425	0.283	0.184	0.150	0.130
19	10	3D	Hx	-	0.851	0.674	0.419	0.296	0.206	0.132	0.099	0.081
20			Sq	-	0.845	0.672	0.418	0.296	0.206	0.132	0.099	0.081
21		2D	Hx	0º	0.857	0.685	0.430	0.311	0.229	0.157	0.125	0.105
22			Sq	0º	0.852	0.681	0.430	0.311	0.228	0.157	0.123	0.104
23			Hx	30º	0.855	0.689	0.428	0.312	0.228	0.157	0.124	0.104
24			Sq	45º	0.862	0.681	0.429	0.313	0.228	0.157	0.125	0.105

Table 5:
Limit stresses in [GPa] for unconstrained load case A and σ_y = 100 MPa.

Nº	Vf [%]	Imperfection	Array	delta	Misalignement angle a
					0.01º	0.1º	1º	2.5º	5º	10º	15º	20º
1	70	3D	Hx	-	4.839	3.624	1.787	1.019	0.607	0.340	0.249	0.202
2			Sq	-	5.162	3.758	1.647	1.008	0.607	0.348	0.259	0.216
3		2D	Hx	0º	4.889	3.652	1.875	1.134	0.686	-	-	-
4			Sq	0º	5.157	3.701	1.722	1.008	0.617	-	-	-
5			Hx	30º	4.880	3.641	1.824	1.068	0.631	-	-	-
6			Sq	45º	5.293	3.973	1.932	1.098	0.676	0.432	-	-
7	50	3D	Hx	-	2.964	2.364	1.324	0.853	0.536	0.312	0.229	0.185
8			Sq	-	2.986	2.292	1.294	0.847	0.537	0.315	0.231	0.188
9		2D	Hx	0º	2.982	2.300	1.346	0.872	0.590	0.372	0.308	-
10			Sq	0º	2.984	2.368	1.325	0.871	0.559	0.356	0.292	-
11			Hx	30º	2.978	2.392	1.369	0.887	0.563	0.365	-	-
12			Sq	45º	3.028	2.435	1.402	0.873	0.580	0.373	0.312	-
13	30	3D	Hx	-	1.928	1.576	0.978	0.681	0.466	0.289	0.215	0.175
14			Sq	-	1.925	1.577	0.977	0.680	0.467	0.290	0.216	0.176
15		2D	Hx	0º	1.943	1.594	1.015	0.717	0.504	0.333	0.272	0.236
16			Sq	0º	1.931	1.570	0.997	0.705	0.496	0.327	0.263	0.227
17			Hx	30º	1.943	1.585	1.006	0.710	0.498	0.330	0.271	0.237
18			Sq	45º	1.929	1.601	1.001	0.719	0.502	0.333	0.275	0.24
19	10	3D	Hx	-	1.134	0.936	0.635	0.480	0.358	0.246	0.192	0.160
20			Sq	-	1.132	0.930	0.637	0.480	0.358	0.246	0.192	0.160
21		2D	Hx	0º	1.138	0.937	0.654	0.508	0.391	0.286	0.234	0.201
22			Sq	0º	1.137	0.944	0.659	0.507	0.391	0.286	0.233	0.199
23			Hx	30º	1.139	0.944	0.654	0.508	0.391	0.287	0.234	0.202
24			Sq	45º	1.138	0.944	0.660	0.508	0.392	0.287	0.235	0.203

Table 6:
Limit stresses in [GPa] for constrained load case B and σ_y = 48.26 MPa.

Nº	Vf [%]	Imperfection	Array	delta	Misalignement angle a
					0.01º	0.1º	1º	2.5º	5º	7.5º	10º	15º	20º
1	70	3D	Hx	-	4.958	3.388	1.303	0.637	0.381	-	-	-	-
2			Sq	-	5.156	3.363	1.145	0.633	-	-	-	-	-
3		2D	Hx	0º	5.001	3.451	1.350	0.741	0.437	-	-	-	-
4			Sq	0º	5.203	3.387	1.257	0.675	-	-	-	-	-
5			Hx	30º	5.021	3.419	1.302	0.679	-	-	-	-	-
6			Sq	45º	4.559	3.242	1.375	0.730	-	-	-	-	-
7	50	3D	Hx	-	3.138	2.301	1.023	0.592	0.351	-	-	-	-
8			Sq	-	3.145	2.267	0.992	0.575	0.353	-	-	-	-
9		2D	Hx	0º	3.171	2.345	1.088	0.630	0.387	-	-	-	-
10			Sq	0º	3.143	2.281	1.040	0.602	0.368	-	-	-	-
11			Hx	30º	3.155	2.322	1.065	0.608	0.367	-	-	-	-
12			Sq	45º	3.224	2.346	1.102	0.625	0.375	-	-	-	-
13	30	3D	Hx	-	2.163	1.590	0.834	0.526	0.338	-	-	-	-
14			Sq	-	2.159	1.590	0.834	0.526	0.339	-	-	-	-
15		2D	Hx	0º	2.178	1.662	0.879	0.553	0.357	-	-	-	-
16			Sq	0º	2.163	1.628	0.850	0.536	0.354	-	-	-	-
17			Hx	30º	2.171	1.659	0.858	0.536	0.351	-	-	-	-
18			Sq	45º	2.174	1.655	0.870	0.549	0.356	-	-	-	-
19	10	3D	Hx	-	1.656	1.288	0.771	0.538	0.380	-	-	-	-
20			Sq	-	1.658	1.288	0.774	0.539	0.381	-	-	-	-
21		2D	Hx	0º	1.666	1.307	0.790	0.538	0.386	-	-	-	-
22			Sq	0º	1.659	1.296	0.782	0.536	0.386	0.317	-	-	-
23			Hx	30º	1.665	1.305	0.787	0.541	0.385	0.317	-	-
24			Sq	45º	1.668	1.310	0.788	0.541	0.386	0.318	-	-	-

Table 7:
Limit stresses in [GPa] for constrained load case B and σ_y = 100 MPa.

Nº	Vf [%]	Imperfection	Array	delta	Misalignement angle a
					0.01º	0.1º	1º	2.5º	5º	7.5º	10º	15º	20º
1	70	3D	Hx	-	6.337	4.719	2.142	1.227	0.723	-	-	-	-
2			Sq	-	6.019	4.390	2.065	1.170	0.726	-	-	-	-
3		2D	Hx	0º	6.392	4.805	2.316	1.351	0.806	-	-	-	-
4			Sq	0º	6.468	4.674	2.153	1.234	0.752	-	-	-	-
5			Hx	30º	6.378	4.765	2.240	1.253	0.738	-	-	-	-
6			Sq	45º	6.220	4.629	2.331	1.341	0.798	-	-	-	-
7	50	3D	Hx	-	3.874	3.014	1.682	1.041	0.652	-	-	-	-
8			Sq	-	3.911	3.079	1.633	1.035	0.654	-	-	-	-
9		2D	Hx	0º	3.933	3.111	1.747	1.110	0.689	-	-	-	-
10			Sq	0º	3.937	3.042	1.687	1.041	0.669	-	-	-	-
11			Hx	30º	3.932	3.125	1.716	1.065	0.671	-	-	-	-
12			Sq	45º	3.995	3.188	1.795	1.111	0.690	-	-	-	-
13	30	3D	Hx	-	2.667	2.168	1.31	0.901	0.615	0.482	-	-	-
14			Sq	-	2.666	2.159	1.32	0.906	0.612	0.482	-	-	-
15		2D	Hx	0º	2.643	2.206	1.34	0.937	0.637	-	-	-	-
16			Sq	0º	2.671	2.187	1.33	0.905	0.624	0.497	-	-	-
17			Hx	30º	2.682	2.159	1.29	0.911	0.624	0.497	-	-	-
18			Sq	45º	2.693	2.194	1.36	0.933	0.629	-	-	-	-
19	10	3D	Hx	-	2.103	1.621	1.144	0.852	0.641	-	-	-	-
20			Sq	-	2.108	1.622	1.145	0.853	0.641	-	-	-	-
21		2D	Hx	0º	2.113	1.547	1.137	0.857	0.651	0.550	-	-	-
22			Sq	0º	2.113	1.547	1.159	0.857	0.649	0.549	-	-	-
23			Hx	30º	2.118	1.587	1.137	0.856	0.650	0.549	-	-	-
24			Sq	45º	2.114	1.588	1.157	0.857	0.650	0.549	-	-	-

3.1 Equilibrum Paths

Equilibrium paths are shown in Figure 7 for selected cases of Groups 4A and 4B (V_f = 70%, 2D imperfection, square fiber array, δ = 0) and σ_y = 48.26 MPa, whereas values of misalignment angles are given in each case. Notice that because only compressive behavior is of interest for micro-buckling, negative values of strain and stress are reported. For small deviation angles, the equilibrium path is fairly linear and reaches a maximum value at a bifurcation load; then the path drops in the post-buckling path. For large deviation angles, the equilibrium path exhibits nonlinearity until a maximum is reached; this is a limit point in the nomenclature of the theory of elastic stability (see, for example, Godoy, 2000); then, the path drops in an unstable behavior.

Figure 7:
Equilibrium paths for selected cases of Group 4 as a function of misalignment angle α. Results for σ_y = 48.26 MPa. (a): Unconstrained load case A; (b): Constrained load case B.

A comparison of load cases A and B shows that the slope is slightly higher in cases B, because the transverse strain has been constrained.

Equilibrium paths for various misalignment angles α, for Groups 4A and 4B and σ_y = 48.26 MPa, are plotted in Figure 8. The scale has been modified with respect to Figure 7 to highlight details of the curves.

Figure 8:
Details of equilibrium paths for selected cases of Group 4 as a function of misalignment angle. Results for σ_y = 48.26 MPa. (a): Unconstrained load case A; (b): Constrained load case B.

For configurations with moderate imperfections α there are local maximum points in the equilibrium paths. But for large values of α the maximum vanishes, with the consequence that the stress increases with a compressive strain without crossing a singularity in the path; micro-buckling does not occur in such cases and the problem is governed by large displacements. This behavior of having a limiting angle α for which the maximum ceases to exist, occurs in most groups studied and is typical of shell buckling problems (Godoy, 2000).

Values of the maximum loads (a limit point in the equilibrium path) in each case are shown in Tables 4 to 7. In groups for which a value is missing means that a maximum in the path was not found, i.e. micro-buckling did not occur. For all cases in Tables 4 to 7 a value of β = 50 was assumed. Such value was selected based on results of experimental tests reported in the literature. Moreover, to illustrate the influence of β, cases with β between 12.5 and 200 were solved for Groups 1, 7, 13, and 19 (3D imperfection and Hx fiber array), load case B, α = 0.01º and σ_y = 48.26 MPa. Results of maximum stresses are shown in Figure 9 and they show a small change in limit stress for high β values.

Figure 9:
Limit stresses as a function of ratio β = Lf/Df. Results for constrained load case B, σ_y = 48.26 MPa, 3D imperfection, Hx array (Groups 1, 7, 13, and 19).

Equilibrium paths for groups 19A and 19B (V_f = 10%, 3D imperfection, hexagonal fiber arrangement Hx, and σ_y = 48.26 MPa) are plotted in Figure 10, for cases with α = 0.01º, 2.5º, and 20º.

Figure 10:
Equilibrium paths for selected cases from Group 11 for different misalignment angles α. Results for σ_y = 48.26 MPa. (a): Unconstrained load case A; (b): Constrained load case B.

In the almost-perfect case, α = 0.01º, a line is shown with the initial slope. There is a change in slope during the loading process of the UC which is caused by plasticity in the matrix. This change in slope occurs at a strain σ₁₁ = -0.012, and is more evident in load case A than in case B. Similar behavior was found for the other cases with V_f = 10%, but the slope change caused by plasticity is not clearly observed in higher volume fractions because the matrix has less incidence in the stiffness of the UC.

3.2 Parametric Studies

Parametric studies have been performed by taking into account the limit stresses in Tables 4 to 7, and results are presented (as in many stability problems) in terms of a knock-down factor η, i.e. the critical stress for a given α divided by the stress for αzero. Because of space restrictions, not all values from Tables 4 to 7 can be plotted.

The knock-down factor η has been plotted as a function of misalignment α for unconstrained load cases A, and constrained cases B, in Figure 11, for 3D imperfections with Sq and Hx fiber arrangement, and several values of misalignment α. The results indicate that the fiber arrangement does not have an incidence on the results. A similar graph for 2D imperfections could also be plotted from values of Tables 4 to 7, and similar results are found.

Figure 11:
Influence of fiber arrangement, considering square (Sq) and hexagonal (Hx). For 3D imperfections knock-down factors η are given for Groups 1, 2, 7, 8, 13, 14, 19, and 20, with σ_y = 48.26 MPa. (a): Unconstrained load cases A; (b): Constrained load cases B.

The incidence of V_f for various misalignment angles is shown in Figure 12, for a given type of imperfection and fiber arrangement (3D and Hx array), and yield stress σ_y = 48.26 MPa. The results show high imperfection-sensitivity to small amplitude imperfections: as an example, for α = 1º, the knock-down factor falls to less than 50% in all cases. This sensitivity increases with increasing V_f . Similar behavior was obtained for 2D imperfections and for Sq configurations and σ_y = 100 MPa.

Figure 12:
Micro-buckling sensitivity. Knock-down factor for Groups 1, 7, 13 and 19 (3D imperfection, Hx array, and σ_y = 48.26 MPa). (a): Unconstrained load cases A; (b): Constrained load cases B.

The influence of the type of imperfection (2D or 3D) has been investigated for groups 1, 3, 7, 9, 13, 15, 19, and 21, considering an hexagonal fiber configuration (Hx) and σ_y = 100 MPa. The knock-down factors are shown in Figure 13, and it is clear that the curves with different types of imperfection become almost identical. The conclusion is that the type of assumed imperfection (either 2D or 3D) does not have a significant influence on the knock-down factor.

Figure 13:
Influence of imperfection type: two-dimensional (2D) and three-dimensional (3D). Knock-down factors for Groups 1, 3, 7, 9, 13, 15, 19 and 21 (Hx array and σ_y = 100 MPa). (a): Unconstrained load cases A; (b): Constrained load cases B.

The influence of angle α (the orientation of the 2D imperfection plane), for Groups 3, 5, 9, 11, 15, 17, 21, and 23, for both load configurations A and B and for the same fiber arrangement (Hx) and σ_y = 100 MPa has been plotted in Figure 14. The results indicate that δ does not have an influence on the limit stress. The same behavior was found for square fiber arrangement (Sq) and σ_y = 48.26 MPa.

Figure 14:
Influence of angle δ (orientation of 2D imperfection plane). Knock-down factors for cases 3, 5, 9, 11, 15, 17, 21 and 23 (2D imperfection, σ_y = 100 MPa, and Hx array) for δ = 0º and 30º, as shown in legend. (a): Unconstrained load cases A; (b): Constrained load cases B.

Cases with higher matrix yield stress produced higher limit stresses. In Figure 15 results for cases 4 and 22 (2D imperfection, Sq array, δ = 0, and σ_y = 48.26 MPa) are shown and the same trend was found in all cases considered in this work. The influence of matrix yield stress on knock-down factor is shown in Figure 16 for cases 4, 16, and 21, with 2D imperfection, Sq fiber arrangement, and δ = 0. The results show an increase in knock-down factor but the problem is still highly sensitive to imperfection.

Figure 15:
Influence of matrix strength on limit stress for cases 4 and 22 (2D imperfection and Sq array). (a): Unconstrained load cases A; (b): Constrained load cases B.

Figure 16:
Influence of matrix strength on knock-down factors for cases 4, 16, and 22 (2D imperfection and Sq array). (a): Unconstrained load cases A; (b): Constrained load cases B.

3.3 Comparison with a Simplified Model

Comparison of the present Finite Element results with a simplified model presented by Barbero (1998) is performed in this section. The analytical equation (eq. 4.93 in Barbero, 2010) is given by

(13)

This analytical equation approximates the limit stress σ(α) of a unidirectional composite under compressive load, without taking into account transverse effects. The composite properties required in this model are the in-plane shear strength F₆ and the in-plane shear modulus G₁₂ . Such results are compared in Figure 17 for groups 1 and 19 (3D configuration, Hx array, and σ_y = 48.26 MPa). For this numerical example and based on the assumed von Mises yield criterion for the matrix parameters, a value F₆ = 27.86 MPa was adopted in concordance with σ_y = 48.26 MPa. The elastic modulus G₁₂ was obtained using the Periodic Microstructure Model (PMM), equation (4.39) in Barbero (2010).

Figure 17:
Comparison with analytical results. Sensitivity curves for Groups 11 and 1 (3D imperfection, σ_y = 48.26 MPa, and Hx array) and analytical equation (Barbero, 2010). (a): Unconstrained load cases A; (b): Constrained load cases B.

The results of the present Finite Element and the approximate analytical model are in very good agreement. Constrained cases A yield lower limit stresses. Similar trends are obtained for all cases considered.

As a reference, the limit stress for both load cases considered has been normalized in Table 8 with respect to bifurcation loads obtained from the model due to Rosen (1965).

Table 8:
Limit stresses normalized with respect to bifurcation loads due to Rosen (1965). Results for α = 0.01º.

	σy = 48.26 MPa		σy= 100 MPa
Vf [%]	Case A	Case B	Case A	Case B
10	0.529	1.029	0.704	1.306
30	0.715	1.045	0.931	1.289
50	0.794	1.082	1.022	1.336
70	0.767	1.027	1.002	1.312

Only results of cases from Group 1, 7, 13, and 19 (3D imperfection and Hx array) for α = 0.01º are considered, and it is found that unconstrained cases A have lower values than constrained cases B. Such discrepancies between simplified model due to Rosen and the present one could be caused by different assumptions in stress states for the fiber and matrix, consideration of fiber misalignment and matrix plasticity, among other reasons. Similar trends were found with Sq array or 2D imperfection.

Finally, Figure 17 and 15 show that there is a loss in fiber reinforcing contribution for large misalignment angles. Also, for such misalignment angles, a change in the matrix yield stress produces limit stresses that are comparable and even higher than those generated by changing V_f . As an example, case 22A in Table 4 (10% of fiber volume fraction, 2D imperfection, σ_y = 48.26 MPa, and Sq array) with α = 5º has a limit stress of 0.228 GPa. If V_f is increased up to 70% (case 4A from Table 4 is recovered), then a value of 0.335 GPa is obtained. If case 22A with α = 5º in Table 4 is considered again and the matrix yield stress is changed from 48.26 MPa to 100 MPa, recovering case 22A from Table 5, a value of 0.391 GPa is obtained. Similar results were found for all cases considered here. Finally, the results seem to show that, for large misalignment angles, the matrix yield strength have a more important role than fiber volume fraction on limit stress.

4 CONCLUSIONS

The micro-buckling of unidirectional fiber-reinforced composites has been investigated in this research by means of a computational micromechanics simulation in which misalignment imperfections were geometrically represented. Results were obtained by means of a Finite Element discretization of the periodic Unit Cell domain, using the general purpose package ABAQUS, and assuming non-linear kinematic and material behavior.

Based on the results, it is possible to conclude that there is a high imperfection-sensitivity in the critical stress of fiber micro-buckling. This is in agreement with earlier results by Drapier et al. (1998) for a unidirectional carbon fiber/epoxy resin material for a fixed imperfection wavelength. A sharp drop is seen to occur for very small angles less than α = 1º, which is consistent with the 2/3 power law identified by Tomblin et al. (1997). For higher values of imperfection amplitude, the asymptotic model of Tomblin et al. (1997) is only an approximation, and more refined values are shown in this work, for example in Tables 4 to 7.

The type of imperfection (either 2D or 3D), the fiber configuration (hexagonal or square), and angle of the plane of 2D imperfection (orientation of 2D imperfection) do not have a great influence on the limit stresses due to micro-buckling. For small values of misalignment, the problem is mainly influenced by fiber volume fraction and matrix yield stress.

For large imperfection amplitude, as given by large misalignment angle, there is a significant loss in the reinforcing effect that is contributed by the fiber. Also, for such imperfection level, increasing the matrix yield stress produces comparable or even higher limit stresses than those produced by fiber volume fraction changes. In other words, micro-buckling seem to be dominated by fiber volume fraction for small misalignment angles, whereas, for large fiber angles it shows as a property dominated by matrix yield stress. A sequel of this conclusion is that increasing fiber volume fraction is not an effective way to increase micro-buckling capacity unless a misalignment control is introduced during the fabrication process.

This work has been restricted to an analysis at the micro level, and no attempt has been made to couple micro and macro levels; this is seen as a topic for further research. The research reported aims to highlight which micromechanics variables play the most important influences on the micro-buckling phenomenon, which is a necessary ingredient before proceeding to multi-scale coupling.

References

• Abaqus v. 6.7, (2009). Dassault Systèmes, Providence, RI, USA.
• Barbero, E. J. (1998). Prediction of compression strength of unidirectional polymer matrix composites. Journal of Composite Materials, 32(5): 483-502.
• Barbero, E. J. (2010). Introduction to composite materials design.CRC Press. Boca Raton FL.
• Budiansky, B., and Fleck, N. A. (1993). Compressive failure of fibre composites. Journal of the Mechanics and Physics of Solids, 41(1): 183-211.
• Car, E., Zalamea, F., Oller, S., Miquel, J., and Oñate, E. (2002). Numerical simulation of fiber reinforced composite materials-two procedures. International Journal of Solids and Structures, 39(7): 1967-1986.
• Drapier, S., Gardin, C., Grandidier, J. C., and Potier-Ferry, M. (1996). Structure effect and microbuckling. Composites Science and Technology, 56(7): 861-867.
• Drapier, S., Grandidier, J. C., and Potier-Ferry, M. (1998). A non-linear numerical approach to the analysis of microbuckling. Composites Science and Technology, 58(5): 785-790.
• Godoy, L. A. (2000). Theory of Elastic Stability: Analysis and Sensitivity. CRC Press. Boca Raton FL.
• Goudarzi, R. H., and Khedmati, M. R. (2015). An experimental investigation of static load capacity of AL-GFRP adhesively bonded single lap and double butt lap joints. Latin American Journal of Solids and Structures, 12(8): 1583-1594.
• Jochum, C., and Grandidier, J. C. (2004). Microbuckling elastic modelling approach of a single carbon fibre embedded in an epoxy matrix. Composites science and technology, 64(16): 2441-2449.
• Jones R.M., Mechanics of Composite Materials, Hemisphere, New York, 1975.
• Kueh, A. B. H. (2013). Buckling of sandwich columns reinforced by triaxial weave fabric composite skin-sheets. International Journal of Mechanical Sciences 66, 45-54.
• Kueh, A. B. H. (2014). Size-influenced mechanical isotropy of singly-plied triaxially woven fabric composites. Composites Part A: Applied Science and Manufacturing 57, 76-87.
• Kueh, A. B. H. and Pellegrino, S. (2007) ABD matrix of single-ply triaxial weave fabric composites. 48th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference, 23-26 April 2007, Honolulu, Hawaii, AIAA-2007-2161.
• Li, S., and Wongsto, A. (2004). Unit cells for micromechanical analyses of particle-reinforced composites. Mechanics of Materials, 36(7): 543-572.
• Maewal, A. (1981). Postbuckling behavior of a periodically laminated medium in compression. International Journal of Solids and Structures, 17(3): 335-344.
• Mallikarachchi, H. M. Y. C. and Pellegrino, S. (2011). Quasi-static folding and deployment of ultrathin composite tape-spring hinges. Journal of Spacecraft and Rockets 48(1), 187-198.
• Mallikarachchi, H. M. Y. C. and Pellegrino, S. (2013). Failure criterion for two-ply plain-weave CFRP laminates. Journal of Composite Materials 47(11), 1357-1375.
• Oller, S., Miquel Canet, J., and Zalamea, F. (2005). Composite material behavior using a homogenization double scale method. Journal of Engineering Mechanics, 131(1): 65-79.
• Rasin, N., Kueh, A. B., Mahat, M. N., and Yassin, A. Y. M. (2016). Stability of triaxially woven fabric composites employing geometrically nonlinear plate model with volume segmentation ABD constitution. Journal of Composite Materials, 50(19) 2719-2735.
• Rosen, B. W. (1965). Fiber composite materials. American Society for Metals, Metals Park, Ohio, 37.
• Sun, C. T., and Tsai, J. L. (2001). Comparison of microbuckling model and kink band model in predicting compressive strength of composites. In Proceedings of the 13th International Conference on Composite Materials, Beijing, China.
• Tomblin, J. S., and Barbero, E. J. (1997). Statistical microbuckling propagation model for compressive strength prediction of fiber-reinforced composites. In Composite Materials: Testing and Design, Thirteenth Volume. ASTM International.
• Tomblin, J. S., Barbero, E. J., and Godoy, L. A. (1997). Imperfection sensitivity of fiber microbuckling in elastic-nonlinear polymer-matrix composites. International Journal of Solids and Structures, 34(13): 1667-1679.
• Van Den Einde, L., Zhao, L., and Seible, F. (2003). Use of FRP composites in civil structural applications. Construction and Building Materials, 17(6): 389-403.
• vanDijk, N. P. (2016). Formulation and implementation of stress-driven and/or strain-driven computational homogenization for finite strain. International Journal for Numerical Methods in Engineering, 107: 1009-1028. doi: 10.1002/nme.5198.
  » https://doi.org/10.1002/nme.5198
• Yurgartis, S. W., and Sternstein, S. S. (1994). Experiments to reveal the role of matrix properties and composite microstructure in longitudinal compression strength. ASTM Special Technical Publication, 1185: 193-193.
• Zahr-Viñuela, J., and Pérez-Castellanos, J. (2011). Elastic constants and isotropy considerations for particulate metal-matrix composites. A multi-particle, cell-based approach. Composites Part A: Applied Science and Manufacturing, 42(5): 521-533.
• Zohdi, T. I., and Wriggers, P. (2008). An Introduction to Computational Micromechanics.Springer

Publication Dates

• Publication in this collection
  Dec 2016

History

• Received
  19 Feb 2016
• Reviewed
  22 Sept 2016
• Accepted
  11 Oct 2016