Modelling of fibre steered plates with coupled thickness variation from overlapping continuous tows

Previous research has hinted on further improvements of the buckling behaviour of variable-stiffness laminates by incorporating overlaps, resulting in a variable thickness profile that is non-linearly coupled to the steering angles. The present study compares two modelling strategies to consider the variable thickness distribution: 1) as manufactured with discrete thicknesses; and 2) smoothed with a continuous thickness distribution. The as-manufactured discrete thickness created by overlapping tows is obtained by means of virtually manufactured laminates. The smeared approximation is much simpler to implement, whereby the local thickness is a non-linear function of the local steering angle. Linear buckling analyses are performed by means of fast semi-analytical models based on the Ritz method using hierarchical polynomials and classical plate formulation. By assuming a smooth manufacturing mould on one side, a one-sided thickness variation is produced, resulting in non-symmetric laminates for which the mid-plane surface is varied accordingly. Modelling guidelines are provided regarding the use of the smeared model in a study covering a wide range of geometries, loading and boundary conditions. With these guidelines, one can apply the smeared thickness technique in semi-analytical models to reach a correlation within ±5% compared to a costly discrete-thickness finite element model.


Introduction
In the aerospace industry, weight reduction has always been an important goal in order to increase efficiency. This has caused a shift in the materials from wood to metal, until lately with composites, such as in the Boeing 787 or Airbus A350. Owing to their specific strength and stiffness, combined with their mechanical tailoring properties, composite materials have been employed progressively more. The latest composite developments involve the emergence of Tow-Steered Laminates (TSL), also called variable-stiffness or Variable-Angle Tow (VAT) laminates, increasing structural capabilities by incorporating a higher degree of tailoring.
Early use of local fibre steering was performed in order to increase open hole strength [47]. Later on, Biggers et al. [3] strategically positioned uni-directional strips in order to improve the buckling capabilities of composite plates, locally varying the stiffness and benefiting from load redistribution. The addition of these strips resulted in Variable Stiffness Laminates (VSL). Building upon this idea of varying the stiffness distribution locally, VAT, also called steered-tow or fibre-steered laminates were created. VAT laminates consist in the spatial variation of the fibre orientations within a ply, as opposed to straight fibres used in conventional Constant Stiffness (CS) laminates. This difference is illustrated in Fig. 1.
Since, the benefits of VAT have been used and researched numerically into increasing mechanical properties, such as stiffness [19], buckling [18,20,29,46] and post-buckling [36,44,45] capabilities for laminates alone, or along the interaction of stiffeners in a bay design [15,33] and linear buckling of VAT filament-wound cylinders [40]. These benefits have been proven in experiments to increase properties in different applications, such as bending [38] and buckling for cylinders [24,25] and plates [27,30].
The production of fibre steering requires automated manufacturing methods. Automated Fibre Placement (AFP) is a consolidated process, largely used in the manufacturing of composite parts. The robotic arms can be configured to bend the tows while following a reference curved path [6]. Novel manufacturing techniques such as Continuous Tow Shearing (CTS) produces curved fibres by means of shearing [23].
During the AFP process, overlaps appear due to the non-parallel path course of variable tow angle, as illustrated in Fig. 2. If these overlaps are to be avoided to keep a more uniform thickness, the distance between two subsequent tows can be increased to obtain gaps instead [4]. Such gaps become resin rich pockets and may become predominant locations for damage and failure [26]. Subsequent cut-andrestart of adjacent tows is another possibility to minimize the presence of overlaps and gaps, at the price of increasing the manufacturing costs [31]. Instead of gaps or overlaps, the CTS method produces a continuously changing thickness profile [23].
When tows are bent in the AFP manufacturing process, fibres are prone to wrinkling and twisting [2,7], as a direct result of compressive strains applied on the inner side of the curved fibres. To prevent this negative effect, it is a common practise to constraint the path curvature, with early assumptions using a minimum radius of 635 mm [28], regardless of the material properties or manufacturing parameters; whereas recent studies clearly state that the maximum curvature depends on deposition rate, temperature and materials [14]. Clancy et al. [14] used laser-assisted AFP of thermoplastic tows and presented processing parameters that allowed a minimum steering radius as low as 400 mm. On the other hand, the CTS process allows a minimum steering radius of 50 mm [23].
The presence of gaps and overlaps influences the ideal response of a VAT laminate. The effects of gaps has been incorporated by the defect layer method, where the stiffness matrix is reduced based on the relative gap percentage in a Finite Element Method (FEM) element [16]. Instead, Mishra et al. [31] adapted this defect layer method by means of a homogenization of the stiffness matrix as a function of the fibre angular distortion. However, despite showing higher strength and stiffness than conventional laminates, gap designs are less stiff and have a lower first-ply failure than the overlap design counterparts [27], even when mass-normalized, in bending [38] and buckling experiments [30]. However, little literature has investigated numerically the effects of overlaps on the buckling behaviour of VAT laminates.
Wu [41] analyzed the effect of a varying thickness profile due to overlaps created by AFP on a cylinder. After replicating the tow courses and overlap locations, as shown in Fig. 3, Wu used this information in a FEM analysis. He defined each element's stiffness with the thickness and layup information of the virtual model closest to the element's centroid and perform a non-linear buckling analysis. This approach showed to be in good comparison with experiments of the manufactured cylinder with overlaps [42].
Upon the variable thickness phenomena, Castro et al. [10] proposed a smeared thickness formulation to approximate the local thickness due to overlaps in AFP-manufactured laminates. The formulation is based on the tow width approximation from Blom et al. [5], whilst applying a constant-volume constraint.
Groh and Weaver [18] investigated the influence of a varying thickness profile in CTS manufactured laminates, modelling the variation both as a 2D plate and as a 3D curved shell, to assess which one represents best the 3D behaviour determined by means of a model using linear brick elements. They concluded that a curved shell formulation represents better the thickness influence, and also saw the possibility of using double curved thickness variation to increase the buckling performance.
Alternatively, Irisarri et al. [22] and Peeters and Abdalla [35] used variable thickness by means of ply drops along curved fibres to retrieve the stacking sequences after a linear buckling optimization using lamination parameters.
The present study focuses on the use of continuous tows in AFP manufacturing allowing the formation of overlaps. Virtuallymanufactured laminates are used to obtain the locations of the overlaps and layer-wise thickness distribution, according to a given fibre orientation and tow width. The resulting information is a discrete thickness distribution that can be used either in a FEM or a semianalytical model. The potential of the smeared thickness approximation by Castro et al. [10] is assessed for representing the influence of overlaps on the buckling behaviour of the AFP laminates herein studied, and compared with the modelling approach using discrete overlaps. Finally, guidelines on the use of the simpler smeared thickness modelling are given for a wide range of design possibilities and manufacturing parameters.

Virtual manufacturing
This section addresses the idealization of variable stiffness by defining the fibre variation, followed by special considerations regarding the overlaps. Finally, the procedure to obtain the discrete thickness distribution of a laminate and the way to incorporate manufacturing constraints are discussed.

Fibre steering
Two coordinate systems are introduced: a global and a manufacturing one; both illustrated in Fig. 4. The manufacturing coordinate is rotated about the global coordinate system through the ply rotation angle γ. The global coordinate system (x; y) is used to represent the laminate properties and associated loading conditions, whereas the manufacturing coordinate system (x 0 ; y 0 ) is used to define the fibre orientation ϕ and the overlap locations on a ply basis.
A one-dimensional fibre variation is defined in the manufacturing coordinate system along the x 0 axis by means of a Lagrange interpolation scheme, as suggested by Wu et al. [46] and shown in Eq. (1), where the manufacturing angle ϕ is shown in Fig. 4. This interpolation scheme allows for complex angle definitions using only a few control points (CP) per ply. In Eq. (1), there are M control points located at fixed positions x 0 k .
The manufactured fibre orientation ϕðx 0 Þ must be translated to the global coordinate angle θ, also defined as positive in the counterclockwise direction and measured about the x-axis, according to Eq.
(2) and illustrated in Fig. 4. The exact method and considerations to rotate a ply are treated in Section 2.2.
A grid in the global coordinate system is created by N s x N s sampling points to keep track of the thickness and stiffness distribution, necessary for the buckling simulations. For each point in this grid, the following parameters are stored for every ply: • Manufacturing angle, ϕ • Global fibre orientation, θ • Fibre steering curvature, ρ • Local thickness (either discrete or smeared), t local

Ply rotation
One can construct the same global fibre orientation with different ply rotation angles, as exemplified by the two cases of Table 1. Two control points (CP) are used, creating a linear variation of the orientation angle ϕ. Note that in both Case 1 and Case 2 the final global fiber orientation is the same, as given by θ ¼ γ þ ϕ of Eq. (2). However, the non-linear thickness variation due to overlaps is only a function of the manufacturing angle ϕ, as discussed in Section 2.3 and 2.4, resulting in different thickness profiles for Case 1 and Case 2. Until now, the literature only assumed uniform thickness at any location, such that the ply rotation angle γ could simply be added afterwards [19,20,26,31,43,46]. By means of this example, it is clear that when overlaps are considered, the rotation angle γ must be judged concurrently with the manufacturing angles ϕ.
The methodology herein suggested, taking the rotation angle into account, is based on the actual production process: a manufacturing area is assumed, according to the manufacturing coordinate system ðx 0 ; y 0 Þ, which encompasses the complete laminate, represented by the black-dotted area in Fig. 5. This manufacturing area can be viewed as the region where tows, such as the one in Fig. 8, are being deposited next to one another by the robot head until fully covering the laminate.
The manufacturing angles ϕ are interpolated in the manufacturing coordinates system according to Eq. (1). However, as the dimension along the x 0 axis varies depending on the rotation angle γ, the CP locations are expressed as percentages rather than fixed positions. In this  way, the CP can always be placed on the edges of the manufacturing area, which are at least corners of the laminate, such that an interpolation is possible over the complete domain. Thereafter, the global laminate mesh is remapped to the manufacturing coordinate system by means of Eq. (3), which can then be used in the interpolation of Eq. (1). The new x 0 still serves as location for the interpolation, however, the y 0 part is also important in order to know the exact overlap location, as there could only be partial overlapping depending on the manufacturing angle, as explained in Section 2.3.
Finally, a positive ply rotation yields a varying manufacturing angle pattern from the top-left corner to the bottom-right corner, and the opposite for a negative rotation. The rotation of a ply can serve two purposes: either increasing the thickness for a given global fibre orientation; to achieve a global orientation not possible due to the constraint on the manufacturing angles, as discussed later; or a combination thereof.
The rotation strategy, and subsequent steps in this methodology, are demonstrated on a virtual manufacturing example for one ply, whose characteristics are given in Table 2. Note that the CP locations are given as percentages along the manufacturing coordinate system ðx 0 ; y 0 Þ and that for each CP there is a corresponding angle. Four CPs are used, leading to a third-order interpolation of the angle ϕ, according to Eq. (1). In Fig. 6, the fibre orientation is plotted in the manufacturing coordinate system with the plate rotated, where the CP values of −30°and −50°are obtained at the top-left and bottom-right corners. When the plate is viewed in its original global coordinate system, the obtained interpolation pattern is the one shown in Fig. 7, where the ply rotation γ has also been applied to obtain the angle distribution θ.

Overlap location
To evaluate the overlap locations, every tow path is determined and plotted [21] using a given opacity level. When the opaque tows superimpose, the distinct opacity levels can be related to discrete thickness values. This procedure is applied for each individual ply, whose separate thickness can then be added together to obtain the total laminate thickness distribution.
By means of the Lagrange interpolation scheme of Eq. (1), the tow path is related to the fibre orientation through Eq. (4) [5] in the manufacturing coordinate system for a given ply. The tow path consists of a sequence of points x 0 ; y 0 , in which x 0 is obtained using fixed increments Δx 0 , and y 0 is obtained using forward Euler integration, where Δy 0 is calculated using Eq. (5). The step size Δx 0 is defined by the length of the encompassed manufacturing domain along the x 0 axis divided by N s À 1. Once the tow path for a single tow is obtained and plotted, the process is repeated by changing the starting point according to the shift direction, until the complete plate is covered. The starting point is shifted by a distance equal to the tow width, assuming that for any ply a manufacturing angle of 0°will always be present, such that no gaps are present. A single tow of the virtual example is shown in Fig. 8 in the manufacturing coordinate system, with the corresponding fibre orientation.
According to the aforementioned shifting strategy, the manufacturing angles ϕ ¼ 0 result in no overlaps or gaps. However, as soon as ϕ -0 the greater effective width w e , given by Eq. (6) [5], will lead to overlapping adjacent tows. The limit case of a single full overlap is obtained for an effective tow width of w e ¼ 2w, covering half of the previous and subsequent tows. This corresponds to jϕj ¼ 60 , according to Eq. (6). Even more overlaps are appearing as jϕj increases. The manufacturing angles and subsequent overlap formations are shown in Fig. 9, with darker regions representing more thickness and overlaps. The exact amount of overlaps is determined by means of pre-established thresholds of the opacity α value. The total local thickness can be obtained by means of Eq. (7), meaning the lowest nominal tow thickness is obtained where no overlaps are present, and increases as an integer multiplication for each given overlap.
As the overlap patterns are periodic along the shifting axis since no variation is present in that direction, only a restrained plotting surface using the full x 0 length could be used in determining the exact overlap locations. Thereafter, the overlaps could be mapped to the complete rotated plate in the manufacturing area, based on their relative coordinate to this repeat surface. For each manufacturing grid point with transformed x 0 and y 0 coordinates, the ply thickness according to Eq. (7) along the fibre orientation obtained by the interpolation scheme of Eq. (2) is stored. The discrete thickness pattern of the virtual manufacturing example is shown in Figs. 10 and 11.

Smeared thickness
Castro et al. [10] proposed a smeared thickness formulation for AFP-manufactured laminates, where the local thickness t local is approximated by Eq. (8). The thickness has a non-linear dependency with the  manufacturing steering angle ϕ defined in Fig. 12, t local and t tow being respectively the local smeared thickness and nominal tow thickness.
This approximation has several limitations, as firstly it is only valid for jϕj < 90 , since manufacturing angles outside this range will yield a negative thickness which is physically not possible. Besides that, angles outside this range would mean tows are pointing in the negative x 0 direction and hence will not cover the remaining part of the ply. Furthermore, as jϕj approaches 90°, the smeared thickness of a nominal tow behaves asymptotically to infinite as shown in Fig. 12. In reality, the case ϕ ¼ 90 would produce a major thickness buildup equal to the number of subsequent tows, which is neither desired nor will it result in infinite thickness. Nevertheless, the smeared approximation for the current study is restricted to jϕj <¼ 60 alongside the considerations discussed in Section 3.1, which corresponds to an increase of 3% in local thickness for a 1°increment of the manufacturing angle ϕ, as illustrated in Fig. 12. The smeared thickness distribution of the virtual manufacturing example shown in Fig. 11 is given in Fig. 13, where the smeared distribution results in the same volume as the discrete distribution.

Manufacturing constraints
The manufacturing process of variable stiffness laminates is significantly constrained by the minimum steering radius, especially for AFP. In order to achieve feasible tow paths, their curvature ρ must be lower than 1.575 1=m, as discussed in Section 1. To evaluate this constraint, the tow path radius can be precisely calculated by means of Eq. (9) [17], which is related to the curvature ρ as ρ ¼ 1=r s .     Applying this expression in the manufacturing coordinate system, the first derivative of the fibre path (dy 0 =dx 0 ) is the tangent of the fibre orientation as shown in Eq. (4) [5], where ϕðx 0 ; y 0 Þ is directly calculated using Eq. (1). The second derivative d 2 y 0 dx 02 becomes: in which the derivative dϕ=dx 0 can be constructed in a closed form using Eq. (1), leading to the expression given in Eq. (11).
The curvature constraint is evaluated for each ply at every grid point in the laminate mesh, remapped to the manufacturing coordinate system. The fibre variation and steering happens in this coordinate system, and is simply rotated when viewed in the global coordinate system: the rotation does not alter the manufacturing constraints. As soon as any mesh location does not satisfy the curvature limit, then     that particular ply can not be manufactured and a different value for the CP location or CP angle has to be determined in order to render the design feasible.

Linear buckling analysis
In this section, the simulation assumptions are discussed, followed by the derivation of the buckling problem for the semi-analytical problem, and finally the FEM implementation.

Laminate stiffness
The VAT laminate mechanical properties will be simulated on a mesoscale by means of the ABD stiffness matrix [37], whose generic formulation is given in Eq. (12). Here, Q ij is the ply stiffness [37], rotated by the local angle θ in the global coordinate system; z defines the position of the ply through-the-thickness. To consider the varying thickness profile, each ply thickness z k À z kÀ1 is given by Eq. (7) and in case of the smeared approach, the ply thickness is given by Eq. (8).
Due to the overlap considerations, the cross section does not remain symmetric, as during manufacturing the plies are added on only one side, assuming that at the other side there is a smooth mould. Therefore, the geometrical mid plane varies throughout the laminate, creating a varying offset. This phenomenon was simulated by Groh et al. [18], where a variable mid plane was used in the curved shell and 2D plate formulations. For the curved shell a varying radius of curvature was used throughout the laminate together with the varying offset, whereas for the 2D plate only the varying offset for the mid plane was used. Both approaches allowed a geometrically accurate representation of the ply thickness distribution over the VAT laminate, but the 2D plate formulation does not consider the coupling between out-of-plane displacements and in-plane strains due to the presence of a shell curvature created by the thickness buildup.
Both 2D and 3D options proved to be in good agreement for angles up to 60°, beyond which only the shell representation renders a good accuracy [18]. The proposed methodology of the present study uses the full ABD matrix in the buckling problem derivation, and limits the manufacturing angle ϕ to 60°, alongside the thickness increase considerations of Section 2.4. Hence, a 2D plate model with correct midplane offset values is selected to represent the behaviour of the VAT laminates with varying thickness.

Semi-analytical approach
The general buckling problem is derived based on the neutral equilibrium criterion of the total potential energy given in Eq. (13). Following the derivation in [9,12], the general form of Eq. (14) can be obtained.
where K is the constitutive stiffness matrix based on the geometry and stiffness properties of the laminate; K g is the geometric stiffness matrix calculated using the pre-buckling in-plane stresses. Their detailed derivation can be found in Castro et al. [9,[11][12][13], and the implementation is publicly available in a Python package [8]. Both K and K g are integrated numerically, given that the constitutive stiffness represented by the ABD matrix spatially changes due to fibre steering. The pre-buckling in-plane stresses are calculated by means of a static analysis. This analysis stems from the first variation of the total potential energy functional: Using the global coordinate system x; y, the in-plane displacements uðx; yÞ; vðx; yÞ and out-of-plane displacement wðx; yÞ can be approximated as: where matrices S u ; S v ; S w contain the shape functions of the approximated displacements, detailed next; and c contain the corresponding Ritz coefficients. Thus, the variation of the strain energy δU can be written as: The work of the external forces δV, considering only a force vector f acting on the boundary @Ω of the two-dimensional plate domain Ω, can be written as: Replacing these expressions in Eq. (15) leads to: which holds true for any δc. Hence, the pre-buckling static equation is obtained: If the in-plane forces at the boundary are given by known functions f x ðx; yÞ; f y ðx; yÞ; f can be calculated as: For displacement-controlled problems, the shape functions based on the Legendre polynomials allow a precise partition of the Ritz coefficient's vector c such that the displacements and rotations of each edge can be independently defined. For instance, setting flag t1 and flag t2 for u and v in Eq. (23). In this case, the partition of f corresponding to the unknown Ritz coefficients can be calculated from the prebuckling static equation, and the remaining system can be solved for the unknown Ritz coefficients. Displacement-controlled problems are not included in the present work and are left as a topic for future research.
In the present work the displacement field is approximated using shape functions based on Legendre hierarchic polynomials, as derived by Rodrigues [34]. These polynomials have the advantage to allow a full control of the boundary conditions for the out-of-plane displacement using the 4 first terms given in Eq. (22), whereas the higherorder terms given in Eq. (24) are used for enriching the inner domain. For the in-plane displacements, the boundary conditions are controlled using the 2 terms of Eq. (23) instead of Eq. (22), assuming that S u;v i¼3 ¼ S u;v i¼4 ¼ 0. For all cases, the inclusion of flags allows (flag = 1) or prohibits (flag = 0) translations (t 1 and t 2 ) and rotations (r 1 and r 2 ) at the plate edges, enabling the method to work with any combination of simply-supported, clamped and free boundary conditions [9]. ðÀ1Þ p ð2i À 2p À 7Þ!! 2 p p!ði À 2p À 1Þ! ξ ðiÀ2pÀ1Þ ð24Þ Table 3 shows how the shape functions of Eqs. (22)-(24) are ordered for each field variable. The order of the approximation polynomial along each direction ξ; η is controlled separately and the final approximation function is shown in Eq. (25). Details regarding the implementation are given in Castro and Donadon [9].

FEM
The four-node shell element S4 available in ABAQUS [39] was adopted for the finite element buckling analysis, using full numerical integration as recommended in [5,43,46]. The element properties are assigned based on the layup information of the manufacturing grid point closest to the element's centroid, as performed in [41].

Test cases
Three different cases are used to investigated the use of the smeared thickness approximation in buckling simulations. All the cases have four control points (CP) for each ply, located respectively at 0%, 33%, 66% and 100%. A generic layup is applied to each case, as shown in Table 4, according to the notation [γhϕ 1 jϕ 2 jϕ 3 jϕ 4 i] suggested by Guerdal et al. [19]. The orthotropic material properties are E xx = 181e3 GPa, E yy = 10.3e3 GPa, ν xy = 0.28, G xy = 7.2e3 GPa and t tow = 0.127 mm [46], with a tow width w tow of 6.35 mm, or 1/4 inch.
Case 1 Case 1, represented in Fig. 14, is an equal bi-directional loading case for a square plate, simply supported on three sides and pinned on the last, with four plies and dimensions of 1500x1500 mm. Case 2 Case 2, shown in Fig. 15, represents a flange, where a plate of higher aspect ratio is simulated, with three edges simply supported and one long edge free. It is only loaded in the x direction, and has a length of 1000 mm and width of 200 mm. The laminate consists of four plies. Case 3 Case 3, represented in Fig. 16, is a multi loading case being applied to 2000 x 1000 mm geometry. It is loaded by a running load in the y direction 30% of the distributed load in the x direction and a shear load applied at the opposite side of the clamp with a total load half that of the x direction. The laminate consists of eight plies.
The thickness distribution of the obtained patterns for both discrete and smeared simulation are compared for all three cases of Section 3.4 and given in Figs. 18-20. There is a good visual agreement in the thickness patterns for all three cases herein investigated. Comparing the values of the maximum thickness t max in Table 5 it can be noted that, for the discrete case, higher thicknesses are achieved due to the locally overlapping tows; whereas in the smeared case these overlaps are averaged out, resulting in a smoother thickness distribution and therefore a lower value for t max . The quantitative accuracy of the obtained thickness pattern is performed by means the same volume requirement [10], here translated to an average thickness requirement. Table 5 presents the relative average thickness error calculated using the discrete and smeared configurations, showing that a maximum error of 0:14% occurred for Case 1.

Linear buckling results
Before comparing the results of the different simulations, the convergence of the models is checked: this is assessed in the FEM model by changing the number of elements, whereas in the semi-analytical model, the convergence is controlled by the amount of integration points and the order of the Legendre polynomials constituting the shape functions. By default, 200x200 virtual manufacturing grid points N s are used to generate the thickness and layup information for the FEM model, augmented with the integration point locations in case of the semi-analytical model, along the material properties presented in Section 3.4.

Smeared thickness approximation
The convergence of the FEM results with smeared thickness is given in Appendix A, and these results are used predominately for verifying the semi-analytical implementation. Fig. 21 shows the convergence of the first eigenvalue for all three cases using the same amount of integration points in each direction (N) equal to the order of the Legendre polynomials (cf. Table 3); with one extra integration point (N+1); and with three extra integration points (N+3). The idea of evaluating over-integrated schemes came   from the possibility that the variable stiffness requires more spatial discretization than the displacement fields being approximated. The results, however, show that using 12 shape functions and an amount of integration points equal to N suffices in obtaining a converged result.

Semi-analytical modelling
For the semi-analytical simulation of the discrete thickness case, Fig. 22 shows again the convergence results for N, N+1 and N+3 inte-    Fig. 22 that there is more oscillations in the convergence of the discrete thickness models, when compared to the convergence of the smoothed thickness model given in Fig. 21. Increasing the amount of integration points for low-order shape functions, makes the result oscillate up to 20%. This is expected since the discrete thickness generates a discontinuous stiffness, which is poorly approximated by quadrature rules original developed for continuous functions [32]. For higher-order shape functions the solution converges, whilst showing smaller oscillations for the same reason. Future studies could attempt different numerical integration schemes, such as presented in Abedian and Düster [1], to evaluate and mitigate the observed convergence issues. Nevertheless, for the present work, the required number of shape functions and integration points for a converged semianalytical simulation can not be predicted with certainty, and therefore the discrete thickness approach is not advised in the context of semi-analytical modelling.

Finite element method
The FE buckling eigenvalue convergence is shown in Fig. 23 for all three cases. It does not suffer from oscillations, in contrast to the semianalytical formulation. Using a 1% convergence criterion, convergence is attained across all three cases with 50 elements in each direction.
The converged simulations are finally compared to a much more refined virtual manufacturing model and simulation, to assess the influence of the thickness grid defining the element's properties. These refined models are run with a 1000x1000 virtual manufacturing grid and 100 elements in each direction, with the results presented in Table 6. This comparison shows little influence of the grid refinement, and that the selected 200 × 200 manufacturing grid is sufficient to provide the required information on thickness and layup for a converged FEM simulation for the discrete thickness case.

Discussion
The FEM results of the discrete thickness converge with 50 elements in each direction, using a convergence criterion of 1%. On the other hand, it is difficult to express parameters for a converged semi-analytical model using discrete thickness, as it does not approximate the discontinuous thickness well, leading to an oscillation of the results even with 16 th -order Legendre polynomials. Despite this oscillation, the global trend from the semi-analytical model approach the FEM results. Finally, it was verified that a 200 × 200 virtual manufacturing grid is enough to convey the required thickness and layup information for the discrete thickness representation to the FEM simulation. Fig. 21. Semi-analytical model convergence of the smeared thickness approach for N, N+1 and N+3 integration points. Fig. 22. Semi-analytical model convergence of the discrete thickness approach for N, N+1 and N+3 integration points; where N is the order of the Legendre polynomials used as shape functions. Fig. 23. FE convergence using the discrete thickness approach.

Validity of the thickness approximation
Linear buckling analyses performed using the smeared thickness and semi-analytical modeling are further compared to the discrete thickness FEM simulation. Table 7 gives an overview of the first eigenvalue λ obtained for each of the three cases using both ways of simulation. To further assess whether differences are case or layup specific, two additional random layups (given in Appendix B) are simulated both ways for every case. These results can be found in Table 8.   Firstly, considering all the data, the smeared modelling has a AE5% relative error range with respect to its discrete thickness counterpart. This is a satisfactory agreement between the thickness approximation provided by the smeared approach and the real thickness distribution of the discrete model. Secondly, the range of relative error is regarded as being intrinsic to a given layup sequence and not linked to the smeared approximation. It is assumed to stem from the influence of the local varying discrete thicknesses as opposed to the smooth smeared approximation, as presented in Table 5. The patterns of the ABD components distribution coincide, as was the case with the total thickness build-up of Section 4. The locally different thickness can be of significant impact, as the B and D matrices depend respectively on the square and cubic of the thickness, alongside the fibre orientation. Coupled to the complex and interlinked relations forming the linear constitutive and geometric stiffness matrix, this altogether is believed to lead to the errors seen in Table 7 and 8. Finally, the eigenmode of the first eigenvalue of each test case is compared for both modelling of the smeared semi-analytical model, and the discrete thickness FEM analysis in Figs. 24-26, showing similar patterns. Therefore, it is concluded that the smeared thickness approximation suggested by  [10] is well suited to represent the behaviour of its discrete thickness VAT laminate counterpart.

Tow width ratio limit
The aforementioned results show that the smeared thickness approach gives a good approximation for the linear buckling response of VAT laminates with discrete overlaps. However, this simplified thickness profile has limited representation capabilities. The approximation is well suited for global responses, while local phenomena such as stress concentration, local strength and failure cannot be predicted, especially at and around overlap regions. These aspects should be evaluated locally with a refined model of the exact physical representation, as the smeared thickness approximation does not take into account any information about tow interactions.
\Furthermore, the influence of the tow width is investigated on the validity of the smeared approximation, as it is expected to progressively become a worse approximation for larger tow widths, since the overlaps are more distinct as they appear less frequently, and have a larger area. To this end, all three cases of Section 3.4 are evaluated again with a discrete thickness FEM analysis for different tows widths. The results are shown in Table 9, and the relative difference with their respective smeared eigenvalue in Table 10. In order to compare the effect on the different geometries, the tow width is normalised with the smallest dimension of the structure: this dimension will remain the smallest dimension in any rotated manufacturing frame, whilst being the most influenced by the larger tow widths, creating a more discrete and less smooth thickness distribution. The difference for the normalized dimensions is displayed in Fig. 27.
The data in Fig. 27 does not show the expected deterioration between the discrete thickness FEM simulations and smeared semianalytical solution. Despite the lack of the apparent influence of the tow width, the authors do not believe the smeared approach to hold valid for the complete simulated parameter range. As stated previously, the fibre orientation is taken at exactly the grid points from the interpolation: with narrow tows, the grid points will always be close to a tow centerline. However, when large tows are present many grid points can be within one tow width, but the orientation at those grid points is then dictated by the centerline orientation and bending effect through the width of the tow due to steering, rather than the interpolation scheme, and would therefore induce errors.
Therefore, the effect of the tow width on the validity of the thickness approximation cannot be fully described by the current framework, since the exact fibre orientation is badly predicted for large tow width ratios. However, by visually comparing the discrete thickness given in Figs. 28-30 with their respective smeared profile in Figs. 18-20, the authors propose that: the maximum tow width per minimum dimension ratio should be 0.15. With this limit, it is guaranteed an accurate representation of the laminate stiffness and of the thickness distribution, when using the simple smeared approach. Beyond this threshold, one starts loosing global similarities and local particularities can no longer be smoothed out. Further investigations should be carried out with a detailed tow model taking the fibre bending into account, whereby sampling points can be assigned along the tow's transverse direction, at a large distance from the tow centreline.
The guidelines herein proposed, along with the good approximation of the smeared approach, can be helpful for designers to quickly assess VAT laminates with overlaps in a practical way for buckling as demonstrated in this paper, or just to obtain the pre-buckling or any general linear static state, as explained in Section 3.2. Moreover, the combination with the semi-analytical modelling will serve as a basis for future work on the design and optimisation of VAT structures with thickness variation coupled with the steering angles to increase buckling performance.

Conclusion
In this paper, Variable Angle Tow (VAT) laminates are first replicated virtually to obtain the overlap locations along the fibre angle distribution. The manufacturing fibre variation is parameterized by means of a Lagrange interpolation scheme in the manufacturing coordinate system, which is considered alongside the complete rotation with respect to the laminate coordinate system. Subsequently, the tow paths are plotted graphically with a given opacity, allowing to be superimposed and retrieve a discrete thickness profile and overlap locations. This discrete thickness profile is approximated with a con- tinuous smeared approach, only dictated by the manufacturing angle, to simulate linear buckling. A procedure is also integrated to verify the manufacturability of each ply of the VAT laminates.
Furthermore, the one sided thickness profile of the VAT laminates is incorporated in the simulation with an offset in the ABD formulation, albeit limited to a manufacturing angle of 60°. The semianalytical formulation is based on the Kirchhoff-Love plate kinematics and taking all ABD terms into account, where the displacement fields are solved by means of the Ritz method, with Legendre polynomial shape functions incorporating different boundary conditions through enrichment functions. However, before the thickness and layup information could be used in the simulations, both discrete and smeared thickness profiles are compared with each other to verify the virtual manufacturing methodology, showing good agreement, with similar locations and shape of thickness build-up.
With the pre-processing information verified, the buckling analysis and thickness approximation validity was investigated. In a first instance, the semi-analytical model was verified against a Finite Element Method (FEM) model for different geometries, loading and boundary conditions. This was done based on the smeared approximation with a random variable stiffness layup, and showed good agreement between both simulation methods. Thereafter, the discrete thickness profile was also modelled with both the FEM and semianalytical modelling, with the first method experiencing convergence issues due to the discontinuous thickness. The FEM method converged, whose results were in satisfactory agreement with the semi-analytical outcome using the smeared approximation, for both the buckling eigenvalue and corresponding eigenmode.
However, using the smeared approach for buckling simulations is not valid for any arbitrary value of tow width, and therefore a guide-    line is proposed such that the maximum tow width should be smaller than 15% of the smallest plate dimension, based on the limits of the visual correlation between the smeared and discrete thickness profiles. With this guideline, the smeared thickness approach can be used in the initial assessment of VAT laminates with overlap configurations, and on the design and optimisation of VAT structures with thickness variation coupled with the steering angles.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Smeared thickness FEM convergence
The smeared thickness approach is also simulated by means of FEM, to verify the correct implementation of the semi-analytical model. Fig. 31 shows the first buckling eigenvalue for an increased number of elements, seeded equally in both x and y direction, with the thickness and layup information constructed on a 200 × 200 point manufacturing grid.
From this data, a 1% convergence is achieved with 15 elements for all cases, regardless of the dimensions of the plate. These 15 elements simulations are finally compared to a much more refined virtual manufacturing model, to assess the influence of the grid to construct the elements. This refined model is created with a 1000x1000 point grid for the thickness and layup information, and 100 elements in each direction. The results of this comparison, given in Table 11, show that the refinement of the FEM and manufacturing model has little influence on the buckling outcome, meaning a 200x200 grid is enough for the virtual manufacturing.
The comparison between the converged first eigenvalue λ of both the semi-analytical and FEM simulation for the smeared approximation are given in Table 12. This shows that the semi-analytical predictions are in good agreement with the FEM results, verifying the semianalytical implementation. Table 13 presents the two additional layups for each of the 3 cases to compare the validity of the smeared semi-analytical buckling simulation with the discrete FEM counterpart.