Prediction of Shape Distortions in Composite Wing Structures

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


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 have been made to analytically predict shape distortions of different composite systems Rennick & Radford, 1997, Fernlund G., 2005and 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 .

1
Where and 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.

2
Where ÄT is the change in temperature, and are the coefficient of thermal expansion in the in-plane and transverse directions respectively. Also, and 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. 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. 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 thermoelastic 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 springin 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.  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 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 and chemical shrinkage depend on the temperature T and degree of cure X as , , , 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 .

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 .

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 Latin American Journal of Solids and Structures, 2018, 15 11 Thematic Section , e69 5/15 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 .

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.

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 0 . According to equation 9 , the T g attained after the first step cure is 117 0 . 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 0 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.  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.

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.

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 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.

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 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. 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.

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.