Prediction of Shape Distortions in Composite Wing Structures

Abstract

Shape distortions and warpage are a major source of problems for composite manufacturers. These distortions are usually accompanied by built up residual stresses. They can deform a component so that it becomes useless. It also has the capability to reduce the strength of the structure. In this paper, the three-dimensional version of the constitutive model originally proposed by Svanberg and Holmberg is employed to predict the warpage of a wing planform. The model takes into account important mechanisms such as thermal expansion, resin shrinkage and frozen-in strains developed during curing cycles. The model was implemented into ABAQUS Finite Element code as a user subroutine UMAT. The macromechanical properties of each composite layer were predicted using a micromechanics based approach, implemented into MATLAB. Results show that wings with cross ply laminates with reducing thickness along the span experienced more warpage than quasi-isotropic laminates. Furthermore, for wings with equal thickness along the span, the results show that the quasi-isotropic laminates experienced more warpage than cross ply laminates. Lastly, the results show that wings with progressively reducing thickness experience twist that is varying from the wing root to the wing tip while wings with a constant thickness experience twist mainly at the centre of the wing

Keywords
Residual stresses; Shape distortions; Warpage; Composite Wing; Twist

1 Introduction

Distortions and residual stresses are inherently present in composite structures. Distortions are undesirable effects of the manufacturing process since they deviate the structure from nominal geometry. These distortions also have the potential to compromise aerodynamic performance or make difficult the assemblage of structural components. In the same vein, residual stresses are the driving mechanisms for induced distortions. These stresses can significantly affect the strength of composite structures (Faria, 2015). In view of the above, several studies in the literature has focused on the understanding and the prediction of shape distortions and residual stresses due to manufacturing processes.

Several factors, divided into intrinsic and extrinsic factors, dictate the level of shape distortion and residual stresses (Albert, 2002). Intrinsic factors include, anisotropic thermal contraction and expansion, resin shrinkage, matrix type which affect resin shrinkage, fiber type and fiber volume fraction. Extrinsic factors include cure schedule, shape and size of the structure, type of mold and tool part interaction among others (Cinar, Ozturk, Ersoy, & Wisnom, 2013). These factors all have various contributions to shape distortions and residual stresses in the structure. For anisotropic thermal contraction and expansion, deformation occurs from difference in temperature between the highest curing temperature and room temperature. These values along the longitudinal direction are lower than that along the out-of-plane transverse direction. This anisotropy leads to process-induced warpage in composite parts with curvatures, during ramp-down from the process temperature (Kaushik & Raghavan, 2010). For resin shrinkage, the crosslinking of the polymer due to the curing chemical reaction taking place leads to irreversible resin shrinkage in thermoset composites. Frozen-in strains are deformations which are locked-into the part when the part undergoes phase change from the rubbery phase to glassy phase.

Manufacturing deformations usually manifest themselves as spring-in in angled structures while they manifest as warpage in flat and singly curved sections. Corner sections usually come out with a smaller angle than the corresponding angle of the mold and this phenomenon is called spring-in as shown in Figure 1. Warpage on the other hand is the distortion or twist of a plate out of shape as shown in Figure 2. In the literature, several attempts have been made to analytically predict shape distortions of different composite systems (Rennick & Radford, 1997), (Fernlund G., 2005) and (Kim & White, 1997)). (Rennick & Radford, 1997) showed that the spring-in angle, Δθ for corner sections can be predicted from geometry and expansion strains by using equation (1).

$\text{Δθ}=\text{θ}\left[\frac{{\text{ε}}_1-{\text{ε}}_3}{1+{\text{ε}}_3}\right]\approx \text{θ}\left({\text{ε}}_1-{\text{ε}}_3\right)$ (1)

Where ${\text{ε}}_1$ and ${\text{ε}}_3$ are the free expansion in the in-plane and through-thickness directions respectively. As shown in Figure 1, Äè is the spring-in angle and è is the angle surrounded by the bend. In this way, the strains due to thermal expansion and cure shrinkage can be accounted for. Equation (1) can be modified to capture the contributions of thermal expansion and cure shrinkage as described in equation (2) below.

$\text{Δθ}=\left(\text{θ}-{180}^0\right)\left[\frac{\left({\text{α}}_1-{\text{α}}_3\right)\text{ΔT}}{1+{\text{α}}_3\text{ΔT}}+\frac{{\text{β}}_1-{\text{β}}_3}{1+{\text{β}}_3}\right]$ (2)

Where ÄT is the change in temperature, ${\text{α}}_1$ and ${\text{α}}_3$ are the coefficient of thermal expansion in the in-plane and transverse directions respectively. Also, ${\text{β}}_1$ and ${\text{β}}_3$ are the coefficient of chemical shrinkage in the in-plane and transverse directions respectively. Essentially, this equation shows that shape distortion is largely driven by the difference between in-plane strains and out-of-plane strains of a structure.

Figure 1
Spring-in of corner sections and curved panels

There exist in the literature several attempts to characterize warpage and spring-in experimentally ((Svanberg & Holmberg, 2001), (Twigg, 2004), (Mahadik & Potter, 2013), (Ali at al, 2017), (Causse et al, 2012) and (Nawab et al, 2017)). Svanberg & Holmberg (2001) presented an experimental data for spring-in of glass fiber epoxy composites. The experiments were performed with angle brackets manufactured by RTM, in a steel mold with accurate temperature control. Different in-mold temperature were used to point out and separate different mechanisms responsible for spring-in. The investigation revealed that three mechanisms were majorly responsible for the spring-in phenomenon observed. The investigation also revealed that cure schedule affects spring-in. The three mechanisms responsible for the shape distortions are: thermal expansion (different in glassy and rubbery state), chemical shrinkage and frozen-in deformations. A simple spring-in model that accounts for these mechanisms was used to illustrate and separate the contribution from each mechanism at different cure schedules, showing that they are all significant. On the other hand, Ali et al (2017) showed that varying the thickness of parts made by increasing the thickness at different sections have different deformation behaviour as compared to uniform thickness parts. It was also discovered that base point has major role in distortion and increasing the thickness at the base changes the mechanics and ability to resist the deformations is increased. Increasing the thickness in flanges also effectively reduces the warpage contribution of the flanges to the overall deformation. Furthermore, Causse et al (2012) showed that the manufacturing quality have a significant impact on the thermo-elastic distortion through different mechanisms depending on the type of defect. Firstly, thickness gradients change locally the fibre volume fraction which in turn modify the coefficients of thermal expansion in the corners of the composite specimens. Secondly, resin rich zones constitute regions of higher CTE (Coefficient of Thermal Expansion) that create stress concentration on one side of the curved area.

In the same vein, some works in the literature focused on using the finite element method to predict the spring-in and warpage behaviours (Fernlund G. O., 2003), (Twigg, 2004), (Dong, 2009)). Fernlund et al in (Fernlund G. O., 2003) used a 6-step method that uses a two-dimensional (2D) special purpose finite element (FE) based process simulation code and a standard 3D structural FE code. The approach avoids the need to develop a full 3D process model, significantly reducing the computational effort yet retaining much of the detail required for accurate analysis. Twigg in his study found a relatively linear increase in maximum part warpage with increasing part length and a decrease in warpage with decreasing part thickness (Twigg, 2004). He also found a strong effect of pressure and tool surface condition was observed in the experimental results. He concluded that a proportionality exists between part warpage and tool CTE. Dong (2009) investigated the use of effective coefficients of thermal expansion in 3D modelling for spring in prediction. He observed that layer-wise modelling for spring in can be laborious and computationally expensive. Using the effective CTE in 3-D finite element analysis, the results showed good agreement with the layer-wise approach. This approach reduced significantly the number of elements and the computation time.

Figure 2
Warpage of Flat panels

Lastly, there exist in the literature different constitutive modelling approaches that are been used in the prediction of shape distortions of composite structures (Baran et al 2017). The models include the viscoelastic model (Zocher, 1997), (Ding et al, 2015), (Nielsen M. W., 2012)), the CHILE - Cure Hardening Instantaneously Linear Elastic model (Bogetti & Gillespie Jr, 1992), (Johnston A. A., 1997)) and the path dependent model (Svanberg & Holmberg, 2004a), (Ersoy et al., 2010)). In the case of the viscoelastic model (Zocher, 1997), the stiffness of the resin depends on the degree of cure in that the resin curing process is divided into three regions. In the first region, the resin was assumed to be fully uncured and assumed to behave like a viscous fluid. The resin modulus is kept constant. In the second region, the resin modulus is considered to be a function of resin’s degree of cure. Also in this region, chemical shrinkage is assumed to occur. In the last region, no further chemical shrinkage occurs and the resin modulus is kept constant. Based on the modulus of the resin, the overall material mechanical properties are calculated. The viscoelastic model has the ability of presenting accurate expression of the mechanical behaviour of resin-dominated composite structures but requires extensive material data for numerical simulations. In the case of the CHILE model ((Johnston A. A., 1997), (Bogetti & Gillespie Jr, 1992)), the fibre is assumed constant during curing. Therefore the emphases is placed on the evolution of resin properties for CHILE law. However, according to (Ding et al 2016), the required run-time for the CHILE model is significantly larger than that for other laws especially in the greater total time increments because each state variable is required to be updated at every time step. In the path dependent model (Svanberg & Holmberg, 2004a), (Ersoy, et al., 2010), a simple mechanical constitutive model that accounted for thermal expansion, chemical shrinkage and frozen-in deformations was presented. The model is a limiting case of linear viscoelasticity that removes the need for rate dependence and replaces it with path dependence on the state variables: strain, degree of cure and temperature. The model, which was derived in incremental form, is applicable for a homogenized curing composite. The model is quite accurate and simple to implement. Additionally, the required run-time for path-dependent law is less than that for viscoelastic law and the CHILE model (Ding et al 2016). This makes the path dependent model quite versatile.

It is common that we make composite parts such as wind turbines, wing structures and conical structures with varying thickness keeping in view its application and design factors like weight reduction, stiffness requirement among others. No study has been reported on process induced residual stress and deformation behaviour of varying thickness composite parts on wing-like composite structures. The objective of this paper therefore is to numerically study the effect of varying thickness on the process induced distortions of wing-like composite structures. The three-dimensional version of the constitutive model originally proposed by Svanberg and Holmberg is employed to predict the shape distortions. A three dimensional approach was chosen because it gives the freedom to model structures with varying thicknesses and profile while still retaining the necessary accuracy desired.

2 Models

2.1 Svanberg model

The model as described by Svanberg in (Svanberg & Holmberg, 2004a) is based on incremental stresses. The constitutive equation is written as

$\text{Δσ}=\{\begin{array}{c}{\text{C}}_{\text{r}}\Delta \left(\text{ε}-{\text{ε}}_{\text{t}}-{\text{ε}}_{\text{c}}\right)-s,T\ge {\text{T}}_{\text{g}}\left(X\right)\\ {\text{C}}_{\text{g}}\Delta \left(\text{ε}-{\text{ε}}_{\text{t}}-{\text{ε}}_{\text{c}}\right),T<{\text{T}}_{\text{g}}\left(\text{X}\right)\end{array}$ (3)

Where ${\text{C}}_{\text{r}}$ is the rubbery modulus tensor and ${\text{C}}_{\text{g}}$ is the glassy modulus tensor. In the most general situation ∆ó contains six stress components: ∆óxx, ∆óyy, ∆ózz, ∆ôyz, ∆ôxz and ∆ôxy. Similar observation holds for $\text{Δε}$, ${\text{Äε}}_{\text{t}}$and ${\text{Äε}}_{\text{c}}$ each one also containing six strain components. The incremental thermal and chemical strains computed in equations (4) and (5) are

${\text{Äε}}_{\text{t}}=\text{αÄT}$ (4)

$Ä{\text{ε}}_{\text{c}}=âÄX.$ (5)

The state variables s are

$\text{s}\left(\text{t}+\text{Δt}\right)=\{\begin{array}{c}\mathrm{0,}T\ge {\text{T}}_{\text{g}}\left(X\right)\\ \text{s}\left(\text{t}\right)+({\text{C}}_{\text{g}}-{\text{C}}_{\text{r}})\Delta \left(\text{ε}-{\text{ε}}_{\text{t}}-{\text{ε}}_{\text{c}}\right),T<{\text{T}}_{\text{g}}\left(\text{X}\right)\end{array}$ (6)

The state variables correspond to stresses that keep track of the loading (thermal, mechanical or chemical) history. In the rubbery state the state variables become zero, meaning that the stress history has been erased. In the glassy state, the s variable represent frozen-in stresses.

The coefficients of thermal expansion ($\text{α}$) and chemical shrinkage ($\text{β}$) depend on the temperature (T) and degree of cure (X) as

$\text{α}=\{\begin{array}{c}{\text{α}}_{\text{l}},X<{\text{X}}_{\text{gel}},T\ge {\text{T}}_{\text{g}}\left(X\right)\\ {\text{α}}_{\text{r}},X\ge {\text{X}}_{\text{gel}},T\ge {\text{T}}_{\text{g}}\left(X\right)\\ {\text{α}}_{\text{g}},T<{\text{T}}_{\text{g}}\left(\text{X}\right)\end{array}$ (7)

$\text{β}=\{\begin{array}{c}{\text{β}}_{\text{l}},X<{\text{X}}_{\text{gel}},T\ge {\text{T}}_{\text{g}}\left(X\right)\\ {\text{β}}_{\text{r}},X\ge {\text{X}}_{\text{gel}},T\ge {\text{T}}_{\text{g}}\left(X\right)\\ {\text{β}}_{\text{g}},T<{\text{T}}_{\text{g}}\left(\text{X}\right)\end{array}$ (8)

where the subscripts 'l', 'r' and 'g' refer respectively to liquid, rubbery and glassy states. Tg is the glass transition temperature and Xgel denotes degree of cure at gelation. The glass transition temperature relates Tg to X by the DiBenedetto equation (Nielsen L. E., 1969) stated in equation (9).

$\frac{{\text{T}}_{\text{g}}-{\text{T}}_{{\text{g}}_0}}{{\text{T}}_{{\text{g}}_{\infty }}-{\text{T}}_{{\text{g}}_0}}=\frac{\text{λX}}{1-\left(1-\text{λ}\right)\text{X}}$ (9)

The mechanical constitutive model has been implemented in ABAQUS as a user subroutine, UMAT (Abaqus, 2011) by using an 8-node hexahedron linear brick element with reduced integration and hourglass control (C3D8R).

3 Validation of model

3.1 Model description

In order to validate the 3d model, the angled bracket simulated by Svanberg & Holmberg, (2004b) was modelled. The geometry, loading and the materials used were the same as that used by by Svanberg & Holmberg, (2004b). The angle bracket was modelled using the 8-node linear brick, reduced integration, hourglass control element (C3D8R). The bracket was modelled with four elements through the thickness with each element representing 1 mm. Also, one element for each 7.5 degree along the curved section. Taking advantage of symmetry, only half of the structure was modelled. The mesh of the bracket is shown in Figure 3a. In Svanberg analysis, 3 different boundary conditions were used. The three boundary conditions used were the contact boundary condition, the free boundary condition and the fully constrained boundary condition. For the purpose of validation, 2 of the boundary conditions have been modelled (the free boundary condition and the fully constrained boundary condition). For the free boundary condition case, the bracket constrained against rigid body motion while it was free to move during the entire cure schedule. In the fully constrained boundary condition case, all nodes were fully constrained to obtain a crude approximation of the in mould constraint. The simulation was done in 4 steps (in mould cure step, demoulding step, post cure step and the final cool down step).

Figure 3
Comparison of Svanberg’s plane strain results with the 3D results

3.2 Results

The results were validated against the results obtained by Svanberg’s plane strain results. Figure 3b shows the deformation experienced after postcure. Figure 4(a) shows the spring in of the angled structure after the first cure step while Figure 4(b) shows the spring in of the angled structure after the second cure step. The results show that the difference between the 3D results and the plane strain results was between 4-6 percent thereby validating the 3D model.

Figure 4
Comparison of Svanberg’s plane strain results with the 3D results

4 Case Study: Composite Wing Application

The case study was done on the upper skin of a NACA 4415 wing of span 1500 mm with a chord of 600 mm as shown in figure below. The study was done on 2 wing planforms each with a quasi-isotropic layup of [90/45/-45/0]_s and a cross ply lay up of [0/90]_s. The study aims to show the effect of shape distortion on the manufacture of wings and wind turbine blades.

4.1 Model Description

Two different wing models have been simulated. Wing A (first wing) features a wing with 3 regions as shown Figure 5 and Figure 6. The wing planform represents a typical wing planform used on small UAVs. The span of the wing is 1500 mm while the chord is 600 mm. The thickness of the wing reduces by half at every 500 mm as denoted by the regions in the figures. As shown in Figure 6, Region A has a thickness of 4 mm while regions B and C have thicknesses of 2 mm and 1 mm respectively. The laminates lay ups shown in Table 1 were made to be symmetric across the regions with material C being the upper symmetry to material B. Additionally, material B is also the upper symmetry to material A. It should also be noted that material C is also a symmetric laminate.

In contrast to Wing A, Wing B (second wing) features a wing without change of thickness throughout the span. The thickness of Wing B is the average thickness of wing A at 2.3 mm. The laminate layup for wing B are [0/90]_4s and [90/45/-45/0]_2s for the cross ply and Quasi isotropic lay-ups respectively.

For the boundary conditions, Points A and B represent positions where the boundary conditions are imposed against rigid body motion. As explained in (Makinde, Rocha, & Donadon, 2017), four steps are used to simulate the cure process. Step 1 details the in-mold cure step. In this step, the maximum temperature which the part experienced was 100⁰. According to equation (9), the T_g attained after the first step cure is 117⁰. Additionally in this step, the part undergoes 3 state changes. Starting out as a liquid, with no modulus, thermal or chemical shrinkage components. Upon gelation, the part is converted to rubbery state. In the rubbery state, the residual stresses gradually build up. However, as the part is fully constrained, no displacement is observed. Afterwards, the part converts into glassy state when the T_g of the part exceeds the cure temperature. In Step 2, the part was demolded. In Step 3, the part was reheated up to 120⁰ and the matrix cure was completed before the part was cooled. It simulated a free standing position for a part that was demoulded and then post cured. In Step 4, the part was finally released and the final deformation was calculated based on nodal displacements.

Figure 5
Wing Planform showing Boundary conditions Points

4.2 Micromechanical model

The elastic constants, thermal and cure shrinkage coefficients were determined from the mechanical properties of their constituent resin and fibers using the Self Consistent Field Micromechanics (SCFM) equations. The SCFM Equations which are listed in ^{Appendix 1} are given in (Johnston, Vaziri, & Poursartip, 2001), and the AS4/8552 carbon fiber properties are taken from (Bogetti & Gillespie Jr, 1992). The properties of HEXCEL 7781-127 have been taken from (Makinde, Rocha, & Donadon, 2017). Table 2 shows the mechanical properties of the crossply laminates while Table 3 shows the mechanical properties of the quasi-isotropic laminates. Table 4 shows the thermal and cure shrinkage properties of both laminate types.

4.3 Meshing and boundary conditions

The mesh and boundary conditions used for the model is shown in Figure 7. The density of the mesh was to ensure that the each element represents a layer of the woven fabric been simulated through the thickness. Two different mechanical boundary conditions have been used in the model. In the first case which is restricted to step 1, the whole part was fully constrained to simulate the in-mould cure phase of curing and to obtain a crude approximation of the constraint from the mould. In the remaining steps, a boundary condition to prevent rigid body motion was employed. The boundary condition to prevent rigid body motion was restricted to the points shown in the figure at the wing root as described in Table 5.

Figure 6
Wing Planform showing various sections highlighting the thickness and material lay-ups

Table 1
Material lay-up of Wing A

Material	Cross ply Lay-up	Quasi-Isotropic Lay-up
Material A	[0/90]2s	[90/45/-45/0]s
Material B	[0/90]4s	[90/45/-45/0]2s
Material C	[0/90]8s	[90/45/-45/0]4s

Figure 7
Finite Element Mesh and Boundary condition points of the wing root section.

Table 2
Mechanical Properties of HEXCEL 7781/HY5052 and AS4/8552 in a Cross Ply laminate

Symbol	HEXCEL 7781/HY5052		AS4/8552		Units
	Glassy	Rubbery	Glassy	Rubbery
E11 = E22	22900	18690	71570	66100	MPa
E33	8480	2310	11970	3100	MPa
í12	0.09	0.0021	0.04	0.0007	-
í13 = í23	0.45	0.85	0.51	0.87	-
G12	2550	27.5	3360	34.8	MPa
G13 = G23	2430	27.4	3160	34.7	MPa

Table 3
Mechanical Properties of HEXCEL 7781/HY5052 and AS4/8552 in a Quasi Isotropic laminate

Symbol	HEXCEL 7781/HY5052		AS4/8552		Units
	Glassy	Rubbery	Glassy	Rubbery
E11 = E22	12100	14560	43790	44000	MPa
E33	8400	2760	11970	3170	MPa
í12	0.52	0.34	0.4	0.33	-
í13 = í23	0.23	0.55	0.31	0.57	-
G12	6500	5480	18840	16530	MPa
G13 = G23	2430	34.7	3160	34.7	MPa

Table 4
Thermal and Cure shrinkage Properties of HEXCEL 7781/HY5052 and AS4/8552 in the 2 laminate types

Symbol	HEXCEL 7781/HY5052		AS4/8552		Units
	Glassy	Rubbery	Glassy	Rubbery
á11 = á22	15.15	5.49	2.72	-1.12	10e-6/deg Cel
á33	6.65	264.8	60.1	223.98	10e-6/deg Cel
â11= â22	-3.60E-03	-7.91E-05	-1.30E-03	-2.34E-05	%
â33	-0.0217	-0.0351	-0.0198	-0.0294	%

Table 5
Boundary condition imposed at different points

POINT A1	POINT B1	POINT B2
u = v = w = 0	v = w = 0	v = 0

4.4 Results

In order to measure the amount of warpage experienced by the wing, the results were quantified in terms of twist of the whole wing and the displacement of the Mean Aerodynamic Centre (MAC). The twist is the rotation of the wing in the xz plane along the span. The twist was evaluated by comparing the displacement of the nodes on the leading and trailing edges with their former position in terms of degrees. The twist angle was measured by comparing the final location of the leading and trailing edge nodes along the span against their original positions.For the MAC displacement, the displacement in the Z direction of the MAC was plotted to give an indication of the departure from pre simulation position.

Figure 8
Displacement plots of Wing A showing twists and Warpage.

Figure 9
Displacement plots of Wing B showing twists and Warpage.

Figure 8 and Figure 9 shows the results of the simulations. Figure 10(a) shows that the twist deformation of the Wing A (wing with ply drops) while Figure 10(b) shows the twist deformation of the Wing B (the wing with equivalent thickness). Figure 10(a) shows that the crossply laminates for both the AS4 and HEXCEL 7781 experience more twist that the quasi-isotropic laminates. There also exist a point of inflexion in the figure where the twist angle changes. This indicates a double twist of the wing profile. For Figure 10(b), the twist for the crossply laminates also have a higher value than that of the quasi-isotropic laminates. However, the trend of the twist deformation for the equivalent wing is different from that of the other wing. The trend exhibits a steep rise in the twist of the wing from the wing root followed by a steep descent towards the wing tip.

Figure 11 shows the z direction displacement of the MAC of Wing A while Figure 12 shows the z direction displacement of the MAC of Wing B. In Figure 11, the displacement of the quasi-isotropic laminates remained relatively constant all through the span of the Wing A while the displacement of the cross ply laminates experienced a constant increase. However at around region B, the gradient of rise of the cross ply laminate was low. This indicates the point of inflexion of the observed in the twist plots. In Figure 12, the displacement plot of the Wing B shows a continuous increase in the displacement of the MAC. The laminate types all experienced relatively the same increase with exception of the carbon quasi-isotropic laminate.

Figure 10
Twist plots of Wing A and Wing B.

4.5 Discussion

The combined picture painted by the twist plot and the MAC displacement plot shows the distortion been experienced on the wings after manufacture. For Wing A, the result in Figure 10(a) shows that the ply drops have a significant effect on the twist of the cross ply. The twist of the wing was rising up the point of inflexion which is at the center of the wing span. After the point of inflexion, the twist continued to rise. The displacement of the MAC for the Wing A for the cross ply layup was also larger than that of the quasi-isotropic layup as shown in Figure 11. The observation shows that cross ply wings with ply drops experience more distortions than quasi-isotropic laminates with ply drops.

Figure 11
MAC Z direction displacement plots of Wing A

For Wing B, the opposite trend was recorded. The twist plot of Figure 10(b) shows that quasi-isotropic plies experienced more twist than the cross ply layup. However, the twist was confined to the center portion of the wing. The magnitude of the twist experienced by the quasi-isotropic plies is comparable to that of the cross ply laminates in wing A. Furthermore, Figure 12 shows that Wing B experienced a continuous rise in the displacement of the MAC. This confirms that the quasi-isotropic laminates experience more distortions on wing B than cross ply laminates on Wing B.

Figure 12
MAC Z direction displacement plots of Wing B

Lastly, the comparison of the results of both Wing A and B shows that the irrespective of the material of the fibre, the quasi-isotropic laminates have lesser twist when in the Wing A configuration and have a higher twist when in Wing B configuration. Additionally, the crossply laminates experience higher twist in the Wing A configuration irrespective of the material used while they experience lesser twist in the Wing B configuration. The results clearly show that shape distortions have a significant effect on the overall wing coefficient of lift and thereby affecting the overall aerodynamic performance. It also shows that the distortions are dominated by the matrix properties.

5 CONCLUSION

In this paper, the 3D version of Svanberg`s model was used to predict shape distortions on a wing planform. The model takes into account important mechanisms such as thermal expansion, resin shrinkage and frozen-in strains developed during curing cycles of polymer composite components. The proposed model accounts for the effects of the residual stresses induced during manufacturing enabling the design of cost efficient, improved performance and optimized composites aero-structures.

The model was implemented into ABAQUS Finite Element code as a user subroutine UMAT. A study case based on composite wing applications is presented. Results show that wings with cross ply laminates with reducing thickness along the span irrespective of the fibre material used, suffered more warpage than quasi-isotropic laminates with reducing thickness along the span. Furthermore, the results show that the quasi-isotropic laminates experienced more warpage than cross ply laminates for wings with equal thickness along the wing span. The results also show that the distortions are dominated by the properties of the matrix used in the part. Finally, the results show that wings with progressively reducing thickness experience twist that is varying from the wing root to the wing tip while wings with a constant thickness experience twist mainly at the centre of the wing.

6 ACKNOWLEDGEMENTS

This work was performed in Instituto Tecnologico de Aeronatica, Brazil and was funded by the Nigerian Air Force and CNPq (process 301053/2016-2).

REFERENCES

• Abaqus, CAE. (2011). User's manual . Abaqus analysis user’s manual.
• Albert, C. a. (2002). Spring-in and warpage of angled composite laminates. Composites Science and Technology, 62(14), 1895 - 1912.
• Ali, M., Nawab, Y., Saouab, A., Anjum, A. S., & Zeeshan, M. (2017). Fabrication induced spring-back in thermosetting woven composite parts with variable thickness. Journal of Industrial Textiles.
• Baran, I., Cinar, K., Ersoy, N., Akkerman, R., & Hattel, J. H. (2017). A review on the mechanical modeling of composite manufacturing processes. Archives of computational methods in engineering, 24(2), 365-395.
• Bogetti, T. A., & Gillespie Jr, J. W. (1992). Process-Induced Stress and Deformation in Thick-Section Thermoset Composite Laminates. Journal of Composite Materials, 26(5), 626-660.
• Causse, P., Ruiz, E., & Trochu, F. (2012). Spring-in behavior of curved composites manufactured by Flexible Injection. Composites Part A: Applied Science and Manufacturing, 1901--1913.
• Chen, H.-J. a. (1996). Three-dimensional effective moduli of symmetric laminates. Journal of Composite Materials, 30 (6), 906 - 917.
• Cinar, K., Ozturk, U. E., Ersoy, N., & Wisnom, M. R. (2013). Modelling manufacturing deformations in corner sections made of composite materials. Journal of Composite Materials, 48(7), 799 - 813.
• Ding, A., Li, S., Sun, J., Wang, J., & Zu, L. (2016). A comparison of process-induced residual stresses and distortions in composite structures with different constitutive laws. Journal of Reinforced Plastics and Composites, 35(10), 807-823.
• Ding, A., Li, S., Wang, J., & Zu, L. (2015). A three-dimensional thermo-viscoelastic analysis of process-induced residual stress in composite laminates. Composite Structures, 129, 60-69.
• Dong, C. (2009). Modeling the dimensional variations of composites using effective coefficients of thermal expansion. Journal of composite materials, 43(22), 2639 - 2652.
• Ersoy, N., Tomasz, G., Kevin, P., Micheal, R. W., David, P., & Graeme, S. (2010). Modelling of the spring-in phenomenon in curved parts made of a thermosetting composite. Composites Part A: Applied Science and Manufacturing, 41(3), 410- 418.
• Faria, A. R. (2015). Prediction of post cure residual stresses and distortions in the Fabrication of Composite Structures. Sao Jose dos Campos: ITA.
• Fernlund, G. (2005). Spring-in of angled sandwich panels. Composites science and technology, 65(2), 317- 323.
• Fernlund, G. O. (2003). Finite element based prediction of process-induced deformation of autoclaved composite structures using 2D process analysis and 3D structural analysis. Journal of Composite Structures, 62, 223 - 234.
• Johnston, A. A. (1997). An integrated model of the development of process-induced deformation in autoclave processing of composite structures. Ann Arbor Canada: University of British Columbia.
• Johnston, A., Vaziri, R., & Poursartip, A. (2001). A plane strain model for process-induced deformation of laminated composite structures. Journal of Composite Materials, 35(16), 1435-1469.
• Kaushik, V., & Raghavan, J. (2010). Experimental study of tool--part interaction during autoclave processing of thermoset polymer composite structures. Composites Part A: Applied Science and Manufacturing, 41(9), 1210--1218.
• Kim, Y. K., & White, S. R. (1997). Viscoelastic analysis of processing-induced residual stresses in thick composite laminates. Mechanics Of Composite Materials And Structures An International Journal, IV(4), 361-387.
• Mahadik, Y., & Potter, K. (2013). Experimental investigation into the thermoelastic spring-in of curved sandwich panels. Composites Part A: Applied Science and Manufacturing, 49, 68-80.
• Makinde, O. M., Rocha, A. F., & Donadon, M. V. (2017). Prediction of Shape Distortions in Angled Composite Structures. 6th International Symposium on Solid Mechanics, (pp. 243-263). Joinville-Santa Catarina.
• Nawab, Y., Sonnenfeld, C., Saouab, A., Agogué, R., & Beauchêne, P. (2017). Characterisation and modelling of thermal expansion coefficient of woven carbon/epoxy composite and its application to the determination of spring-in. Journal of Composite Materials, 51(11), 1527-1538.
• Nielsen, L. E. (1969). Cross-linking-effects on physical properties of polymers. Journal of Macromolecular Science, Part C, 3(1), 69-103.
• Nielsen, M. W. (2012). Predictions of process induced shape distortions and residual stresses in large fibre reinforced composite laminates. Lyngby, Denmark: Ph. D. thesis, Technical University of Denmark.
• Rennick, T., & Radford, D. (1997). Determination of manufacturing distortion in laminated composite. Proceedings of the 11th International Conference on Composite Materials and Composite Process, (p. 302). Microstruct Woodhead Publishing.
• Svanberg, J. M., & Holmberg, J. A. (2004a). Prediction of shape distortions. Part I. FE-implementation of a path dependent constitutive model. Composites part A: applied science and manufacturing, 35(6), 711 - 721.
• Svanberg, J. M., & Holmberg, J. A. (2004b). Prediction of shape distortions. Part II. Experimental validation and analysis of boundary condition. Composites Part A: Applied Science and Manufacturin, 35 (6), 723 - 734.
• Svanberg, J. M., & Holmberg, J. A. (2001). An experimental investigation on mechanisms for manufacturing induced shape distortions in homogeneous and balanced laminates. Composites Part A: Applied Science and Manufacturing, 32(6), 827--838.
• Twigg, G. a. (2004). Tool--part interaction in composites processing. Part I: experimental investigation and analytical model. Composites Part A: Applied Science and Manufacturing, 35(1), 121 - 133.
• Zocher, M. a. (1997). A three-dimensional finite element formulation for thermoviscoelastic orthotropic media. International Journal for Numerical Methods in Engineering, 40(12), 2267 - 2288.

• Available Online: April 18, 2018

7 Appendix 1

7.1 Calculation of Ply Engineering Constants

Using the self-consistent field model equations given by Bogetti in follows (Bogetti & Gillespie Jr, 1992), the ply mechanical properties are calculated as shown below. The inputs are the transversely isotropic mechanical properties of the fibres (E11f, E33f, G13f, í_12f, and í_23f), the properties of the isotropic resin (Er, ír) and the fibre volume fraction, Vf:

In plane Moduli

${\text{E}}_{11}={\text{E}}_{11\text{f}}{\text{V}}_{\text{f}}+{\text{E}}_{\text{r}}\left(1-{\text{V}}_{\text{f}}\right)+\left[\frac{4\left({\text{ν}}_{\text{r}}-{\text{ν}}_{12\text{f}}^2\right){\text{k}}_{\text{f}}{\text{k}}_{\text{r}}{\text{G}}_{\text{r}}\left(1-{\text{V}}_{\text{f}}\right){\text{V}}_{\text{f}}}{\left({\text{k}}_{\text{f}}+{\text{G}}_{\text{r}}\right){\text{k}}_{\text{r}}+\left({\text{k}}_{\text{f}}-{\text{k}}_{\text{r}}\right){\text{G}}_{\text{r}}{\text{V}}_{\text{f}}}\right]$ (10)

The transverse moduli:

${\text{E}}_{22}={\text{E}}_{33}=\frac{1}{\left(\frac{1}{4{\text{k}}_{\text{T}}}\right)+\left(\frac{1}{4{\text{G}}_{23}}\right)+({\text{ν}}_{12}^2/{\text{E}}_{11})}$ (11)

In plane Shear Moduli:

${\text{G}}_{12}={\text{G}}_{13}={\text{G}}_{\text{r}}\left[\frac{\left({\text{G}}_{13\text{f}}+{\text{G}}_{\text{r}}\right)+\left({\text{G}}_{13\text{f}}-{\text{G}}_{\text{r}}\right){\text{V}}_{\text{f}}}{\left({\text{G}}_{13\text{f}}+{\text{G}}_{\text{r}}\right)-\left({\text{G}}_{13\text{f}}-{\text{G}}_{\text{r}}\right){\text{V}}_{\text{f}}}\right]$ (12)

Transverse Shear Moduli:

${\text{G}}_{23}=\frac{{\text{G}}_{\text{r}}\left[{\text{k}}_{\text{r}}\left({\text{G}}_{\text{r}}+{\text{G}}_{23\text{f}}\right)+2{\text{G}}_{23\text{f}}{\text{G}}_{\text{r}}+{\text{k}}_{\text{r}}\left({\text{G}}_{23\text{f}}-{\text{G}}_{\text{r}}\right){\text{V}}_{\text{f}}\right]}{{\text{k}}_{\text{r}}\left({\text{G}}_{\text{r}}+{\text{G}}_{23\text{f}}\right)+2{\text{G}}_{23\text{f}}{\text{G}}_{\text{r}}-({\text{k}}_{\text{r}}+2{\text{G}}_{\text{r}})\left({\text{G}}_{23\text{f}}-{\text{G}}_{\text{r}}\right){\text{V}}_{\text{f}}}$ (13)

Where

${\text{G}}_{\text{r}}=\frac{{\text{E}}_{\text{r}}}{2\left(1+{\text{ν}}_{\text{r}}\right)}$ (14)

${\text{G}}_{23\text{f}}=\frac{{\text{E}}_{33\text{f}}}{2\left(1+{\text{ν}}_{\text{r}}\right)}$ (15)

The Major poisson’s ratio:

${\text{ν}}_{12}={\text{ν}}_{13}={\text{ν}}_{12\text{f}}{\text{V}}_{\text{f}}+{\text{ν}}_{\text{r}}\left(1-{\text{V}}_{\text{f}}\right)+\left[\frac{{\text{ν}}_{\text{r}}-{\text{ν}}_{12\text{f}})\left({\text{k}}_{\text{r}}-{\text{k}}_{\text{f}}\right){\text{G}}_{\text{r}}\left(1-{\text{V}}_{\text{f}}\right){\text{V}}_{\text{f}}}{({\text{k}}_{\text{f}}+{\text{G}}_{\text{r}}){\text{k}}_{\text{r}}+\left({\text{k}}_{\text{f}}-{\text{k}}_{\text{r}}\right){\text{G}}_{\text{r}}{\text{V}}_{\text{f}}}\right]$ (16)

The transverse poisson’s ratio:

${\text{ν}}_{23}=\frac{2{\text{E}}_{11}{\text{k}}_{\text{T}}-{\text{E}}_{22}{\text{E}}_{11}-4{\text{ν}}_{13}^2{\text{k}}_{\text{T}}{\text{E}}_{22}}{2{\text{E}}_{11}{\text{k}}_{\text{T}}}$ (17)

k_r is the plain strain bulk modulus for the isotropic resin defined by:

${\text{k}}_{\text{r}}=\frac{\text{E}}{2\left(1-\text{ν}-2{\text{ν}}^2\right)}$ (18)

k_f is the plain strain bulk modulus for the fibre is defined by:

${\text{k}}_{\text{f}}=\frac{{\text{E}}_{\text{f}}}{2\left(1-{\text{ν}}_{\text{f}}-2{\text{ν}}_{\text{f}}^2\right)}$ (19)

And k_T is the effective plane strain bulk modulus of the composite defined by

${\text{k}}_{\text{T}}=\frac{\left({\text{k}}_{\text{f}}+{\text{G}}_{\text{r}}\right){\text{k}}_{\text{r}}+\left({\text{k}}_{\text{f}}-{\text{k}}_{\text{r}}\right){\text{G}}_{\text{r}}{\text{V}}_{\text{f}}}{\left({\text{k}}_{\text{f}}+{\text{G}}_{\text{r}}\right)-\left({\text{k}}_{\text{f}}-{\text{k}}_{\text{r}}\right){\text{V}}_{\text{f}}}$ (20)

7.2 Calculation of Ply thermal and Cure Shrinkage coefficients

The coefficients below are used to calculate the thermal or cure shrinkage strains in the lamina. The in plane coefficient:

${\text{α}}_1=\frac{{\text{α}}_{1\text{f}}{\text{E}}_{11\text{f}}{\text{V}}_{\text{f}}+{\text{α}}_{\text{r}}{\text{E}}_{\text{r}}\left(1-{\text{V}}_{\text{f}}\right)}{{\text{E}}_{11\text{f}}{\text{V}}_{\text{f}}+{\text{E}}_{\text{r}}\left(1-{\text{V}}_{\text{f}}\right)}$ (21)

The transverse coefficient:

${\text{α}}_2={\text{α}}_3=\left({\text{α}}_{2\text{f}}+{\text{α}}_{1\text{f}}{\text{ν}}_{13\text{f}}\right){\text{V}}_{\text{f}}+\left({\text{α}}_{\text{r}}+{\text{ν}}_{\text{r}}{\text{α}}_{\text{r}}\right)\left(1-{\text{V}}_{\text{f}}\right)-\left[{\text{ν}}_{13\text{f}}{\text{V}}_{\text{f}}+{\text{ν}}_{\text{r}}\left(1-{\text{V}}_{\text{f}}\right)\right]{\text{α}}_1$ (22)

7.3 3D Model

The 3d laminate model used is as defined in (Chen, 1996). Using the equations below, the 3D stiffness matrix of the whole laminate could be derived based on the stiffness matrix of the contributing laminas.

${\text{C}}_{33}={\left({\displaystyle \sum }_{\text{k}=1}^{\text{N}}\frac{{\text{t}}^{\text{k}}}{{\text{C}}_{33}^{\text{k}}}\right)}^{-1}$ (23)

${\text{C}}_{\text{i}3}={\text{C}}_{33}{\displaystyle \sum }_{\text{k}=1}^{\text{N}}{\text{t}}^{\text{k}}\frac{{\text{C}}_{\text{i}3}^{\text{k}}}{{\text{C}}_{33}^{\text{k}}}\text{i}=\mathrm{1,2,6}$ (24)

${\text{C}}_{\text{ij}}={\displaystyle \sum }_{\text{k}=1}^{\text{N}}{\text{t}}^{\text{k}}{\text{C}}_{\text{ij}}^{\text{k}}+\frac{{\text{C}}_{\text{i}3}{\text{C}}_{\text{j}3}}{{\text{C}}_{33}}-{\displaystyle \sum }_{\text{k}=1}^{\text{N}}{\text{t}}^{\text{k}}\frac{{\text{C}}_{\text{i}3}^{\text{k}}{\text{C}}_{\text{j}3}^{\text{k}}}{{\text{C}}_{33}^{\text{k}}}\text{i},\text{j}=\mathrm{1,2,6}$ (25)

${\text{S}}_{\text{ij}}={\displaystyle \sum }_{\text{k}=1}^{\text{N}}{\text{t}}^{\text{k}}{\text{S}}_{\text{ij}}^{\text{k}}\text{i},\text{j}=\mathrm{4,5}$ (26)

Then the interlaminar shear stiffness C₄₄,C₅₅ and C₄₅ are

$\left[\begin{array}{cc}{\text{C}}_{44}& {\text{C}}_{45}\\ {\text{C}}_{45}& {\text{C}}_{55}\end{array}\right]={\left[\begin{array}{cc}{\text{S}}_{44}& {\text{S}}_{45}\\ {\text{S}}_{45}& {\text{S}}_{55}\end{array}\right]}^{-1}$ (27)

Where h is the total thickness of the laminate, and z is the vertical distance from the midplane, the interlaminar shear compliance S₄₄, S₄₅ and S₅₅ can be obtained by:

${\text{S}}_{\text{ij}}={\displaystyle \sum }_{\text{k}=1}^{\text{N}}\left\{\frac{3}{\text{h}}\left({\text{z}}_{\text{k}}-{\text{z}}_{\text{k}-1}\right)-\frac{4}{{\text{h}}^3}\left({\text{z}}_{\text{k}}^3-{\text{z}}_{\text{k}-1}^3\right)\right\}{\text{S}}_{\text{ij}}^{\text{k}}$ (28)

Publication Dates

• Publication in this collection
  2018

History

• Received
  27 July 2017
• Reviewed
  08 Feb 2018
• Accepted
  02 Apr 2018