Characterization of PLA and design of a 3D printed wing

(1)

DEGREE PROJECT IN AEROSPACE ENGINEERING, SECOND CYCLE, 30 CREDITS

STOCKHOLM, SWEDEN 2018

Characterization of PLA and design of a 3D printed wing

JUAN JOSE SAAMEÑO PEREZ

SUPERVISED BY :

PROF ESSOR ERDAL KAYACAN &

P H D. CANDIDATE YUNUS GOVDELI

(2)

Abstract

This report describes the structural design of a wing for a Vertical Take Off and Landing drone, in which all the structure will be built by fused deposition modeling of polylactic acid (PLA). To perform this design, the material used is first characterized in different orientations using tensile stress tests, Image Correlation and MATLAB. These properties are then input in a MATLAB program specially developed for this project to obtain the optimum skin and spar thickness in the wing for certain flight conditions. Results are finally verified with a 3D model in CAD and scaled wings in bending tests.

(3)

1 Introduction

In this chapter it is first briefly described the additive manufacturing process, its background and the arguments for using this technique to build drones. Finally, the objectives are exposed along with their possible delimitations.

1.1 Background

Fused deposition modeling (FDM) is a manufacturing technique first developed by S. Scott Crump in the late 1980s and first commercialized by Stratasys in 1990 [1, p.124]. However, it became extraordinary popular after 2009, when the patent US 5121329 expired.

This method consists on processing an STL file (STereoLithography file format) mathematically slicing and orienting the model for the building process. Then, small flattened strings of molten material are extruded and deposited on a flat surface, which can be either the tray of the printer or another layer of printed material. As soon as the material leaves the nozzle where it is extruded, it is solified and bonded to the previous printed layer. This results in a piece fully built in thermoplastic (polylactic acid in this case) bottom-up.

1.2 Motivation

The use of Unmanned Aerial Vehicles (UAV) has become exponentially popular during the last years in military, civilian and recreational applications. This has led to investigate and optimize their production. Conventional UAV manufacturing methods (such as injection molding, machining and drilling) result in ex- cessive usage of raw materials and labor force. An alternative to the traditional methods is to employ composite materials, such as carbon fiber composites[2]

because of their high strength to weight ratio. However, their cost is high, they need multiple processing steps and there are some design details which are hard to be executed by this method.

The reasoning behind the choice of additive manufacturing methods to build the drone is as follows. First, it is fast compared to traditional manufacturing methods; complex structures can be uploaded from a CAD model and printed in a few hours enabling a rapid verification and development of design ideas.

Secondly, additive manufacturing machines complete a build in one step with no interaction from the machine operator during the build phase, reducing the dependence on different manufacturing processes, giving the designer greater control over the final product.

A first prototype has already been built for this project with PLA by additive manufacturing. Although the outer structure has already been deeply studied in [3], it is convenient to study new designs for the inner structure of the wing, which could make the UAV lighter. Decreasing the weight of the UAV would suppose less power consumption for batteries, less material used, less cost, easier transportation, etc. A new structure with a new material is studied in this report, and for this, mechanical properties of the new material will be determined to later on be used to dimension the new structure.

(5)

1.3 Objectives

The objectives of this report are:

1. Characterization of the material used.

To obtain the mechanical properties (Elasticity and shear moduli, and Poisson ratios) and study the possible printing orientations.

2. Calculation of skin thickness and weight of the UAV.

3. Verify the calculated model.

(a) Stress on the surface of a 3D model in which loads are simulated on it in a CAD program.

(b) Bending of the wings, applying a load (acting as lift) on 3D printed scaled wings.

1.4 Associated problems and challenges

Below, some difficulties which may be faced while completing the objectives are exposed:

Polylactic acid (PLA) does not have the same mechanical properties in all directions (not isotropic).

Performing approppriate type of tests to obtain the most accurate characteristics.

Avoid the use of gauges, since many tests will have to be performed.

Calculating the loads on the UAV.

Complexity of the model may motivate the use of experimental solu- tions.

Determining the best orientation depending on the loads.

1.5 Delimitation

Facilities and machines to measure properties. Unfortunately, the group does not count with unlimited equipment. Therefore, this will lead to look for a solution in accordance with the school means.

Margin of error during characterization. Although it is inevitable, this margin must be as low as possible, since every mistake made in this phase would be carried along the rest of the project.

There is a minimum printable thickness that the printer can extrude.

(6)

There is also a maximum weight that the UAV can reach, and this is determined, between other factors, the thickness of spars and skin.

Besides, the objective is obtaining a value lower than the current one.

3. Limited time, since all the laboratory work must be finished within 5 months.

(7)

2 Relevant theory and method

In this section, the theory to characterize PLA will be described, along with the steps followed to calculate the weight and thickness of the UAV.

2.1 Characterization of PLA

This is a key process to get accurate and sensible results in further steps. Be- sides, the PLA provider (Shenzen Esun Industrial Co., Ltd) stopped producing transparent PLA in the end of 2017, the material in which the first prototype was built, so it was necessary to purchase a new type of material with unknown properties.

It has been demonstrated that the part orientation during the printing plays an important role in the surface finish, dimensional accuracy, cost and mechanical behaviour. Mechanical properties of AM parts depend on imperfections and manufacturing irregularities in the material, since it is deposited in a way which relies on its bonding to previous layers and filling gaps in between. The closest approach to estimate these properties is the orthotropic model[4].

2.1.1 Orthotropic model

To determine the mechanical properties of an orthotropic material it is necessary to find three Young’s moduli (Ei), three Poisson’s ratios (νij) and three shear moduli (Gij).

Figure 1: Strips from an orthotropic material subjected to uniaxial tension.[5]

According to Hooke’s law, the relationship between stress and small strains is linearly proportional. Considering an orthotropic material as in Figure 1, cutting strips in each of the three directions and subjecting these to a uniform stress in that given material direction, the strains are then:

ε1= σ1

E₁, ε2= −ν12ε1= −ν12

σ1

E₁, and ε3= −ν13ε1= −ν13

σ1

E₁ Similarly in 2− and 3− directions

(8)

ε₂= σ2

E2

, ε₁= −ν₂₁ε₂= −ν₂₁σ2

E2

, and ε₃= −ν₂₃ε₁= −ν₂₃σ2

E2

ε₃= σ₃ E3

, ε₁= −ν₃₁ε₃= −ν₃₁σ₃ E3

, and ε₂= −ν₃₂ε₃= −ν₃₂σ₃ E3

(rule-of-thumb for definition of Poisson ratio: first index refers to action and second index to reaction, i.e., stress in direcion i and strain in direction j is given by ν_ij = −ε_j/ε_i). The total strain in the 1-direction is thus by summation

ε1= σ₁ E1

−ν₂₁σ₂ E2

−ν₃₁σ₃ E3

or written out in full for all components





 ε1

ε2

ε3

γ1

γ₂ γ₃







=







1/E1 −ν21/E2 −ν31/E3 0 0 0

−ν12/E1 1/E2 −ν32/E3 0 0 0

−ν13/E1 −ν23/E2 1/E3 0 0 0

0 0 0 1/G23 0 0

0 0 0 0 1/G₃₁ 0

0 0 0 0 0 1/G₁₂











 σ1

σ2

σ3

τ23

τ₃₁ τ₁₂





 (1) Five specimens for each of the six different orientations have been tested to calculate the nine engineering constants following the method used in Me- chanical property characterization and simulation of fused deposition modeling Polycarbonate parts.[4]. They can be obtained from tensile strength tests as:

E1= ∆σ1

∆ε₁ (2)

ν12=ε2

ε1

(3) The in-plane shear modulus can be calculated from the test of a 45-oriented unidirectional specimen, according to the following equation:

G12= E45°/2 · (1 + ν45°) (4) where 1 is the pulling direction and 2 is the perpendicular direction.

2.1.2 Characterization

To obtain the nine independent constants of the stiffness matrix, thirty samples were tested in six different orientations as shown in Figure 3, five samples corresponding to each orientation. As there are not standard tests for Additive Man- ufacturing parts, samples have been built and tested according to ASTMD638:

Standard Test Method for Tensile Properties of Plastics, because it was preferred instead of ISO for most of the authors studying the mechanical behaviour of AM parts [4]. Type IV coupons were designed with SolidEdge®2017 and converted into STL file to be printed with conversion tolerance 0.001 mm and surface plane angle 10°. Since Type IV samples are smaller than the other types, as shown in

(9)

Figure 2: Specimen dimensions for different types of coupons from ASTM 638.[6]

Figure 2, building time is shorter, and they are more likely to break in the field of view (FOV) where the camera is focused to measure the strain. Otherwise, the type of coupons should not influence the result of the tests [7].

Sets of five coupons were designed with Cubicreator to be built afterwards with a Cubicon Single Plus High-Perfomance 3D Printer and the parameters shown in Table 1. Therefore, six batches (one for each orientation shown in Figure 3) with five coupons in each of them were printed, and they were ready to be tested as the printing process was finished, without any additional operation in most cases. For coupons in direction 5 and 6, support structures were placed to facilitate the printing process, this material had to be removed. Figure 4 shows an example of how the samples were placed in Cubicreator, which is at the same time, how these samples were produced in the printer.

Figure 3: Orientation for tensile test samples.[4]

(10)

Figure 4: Preparation for printing the batch of coupons in orientation 3.

Filament 1.75 mm eSun White PLA Flow 100%

Temperature

Extruder Temp. (160∼260) 210ºC Bed Temp. (40∼120) 65ºC Chamber Temp. (30∼55) 40ºC

Quality

Layer height (0.1∼0.5) 0.2 mm Wall thickness (0.4∼125) 0.4 mm Bottom thickness (0.2∼0.5) 0.2 mm

Sink object (0∼) 0

Support

Fill rate (0∼100) 20%

Distance XY 0.5 mm

Distance Z 0.12 mm

Angle 45º

Infill

Rate (0∼100) 100

Top count (0∼n) 0 ea

Bottom count (0∼n) 0 ea Infill overlap (0∼100) 15 %

Speed

Support speed 200 mm/s Travel speed 200 mm/s Inner wall speed 200 mm/s Outer wall speed 100 mm/s Bottom layer speed 30 mm/s

Infill speed 200 mm/s Retraction speed 40 mm/s

Cool

Fan on height 0.5 mm

Fan min. speed 100 %

Fan max. speed 100 %

Minimal layer time 15 sec Table 1: Printing parameters used for the project

(11)

These settings are set the same as used in the previous drone, because they were proved to be adequate[3] for the drone and one of the main points is to compare the current structure with the one studied in this report, so it would be convenient to have as few changes as possible is this aspect. The infill is 100%, which means that the part is completely filled with fully dense raster toolpaths. The extruded angle filaments in the XY plane are +45°/−45°

alternating in each layer, which is the default movement of the printer because it is also demonstrated that this filling style bears better combined loads [4].

Support material was needed to properly print test samples 5 and 6, and to a lesser extent, sample 2. However, the density in this case is 20%, since it will be removed after built and will not affect the part to test. The contour, finished and rest of parameters have been selected to emulate the printing of the UAV.

The tensile tests were performed using a Shimadzu AGS-X 50 kN Universal Testing Machine (UTM) equipped with a high speed video camera to measure the strain in the gauge area (Figure 5). ASTM standards[6] states to use a testing speed “which gives rupture within 0.5 to 5-min testing time”, so the speed was set to 1 mm/min for all samples excepting the ones oriented in direction 3 (Figure 3), whose test speed was set to 0.5 mm/min for consistency. Samples were breaking within the first twenty seconds, preventing the collection of the same amount of data as in the rest. Before placing the samples in the machine, a random speckle pattern was sprayed on all of them, as seen in Figure 6 so the elongation could be recognized afterwards by a Digital Image Correlation (DIC) software called GOM Image Correlate 2017.

Figure 5: Shimadzu AGS-X 50 kN UTM.

Figure 6: Close-up view of a sprayed coupon.

The sprayed surface helps the software to distinguish the relative displacement of the pixels in a designated area, which corresponds in this case, only to the narrow section of the coupon as seen in Figure 7. Axial and longitudinal elongation was collected in this area as the average of the relative displacement of every point.

It was also necessary to decrease the number of frames per second (fps) of the video from ∼ 24 to 10 fps since the stress data from the UTM is collected ten times per second. This conversion allows to relate the stress data set with its corresponding strain data set to operate with them.

(12)

Figure 7: View of the coupon with speckle pattern seen in the DIC software with strain vs. time in the diagram.

2.1.3 Program ”E Calculator.m”

More than 60 samples were tested (the tests had to be repeated once due to some defect in the first batch) to properly characterize this type of PLA. Therefore, it was convenient to automize as much as possible the process to obtain these properties. A program (see appendix 6.1) was developed in Matlab together with Enea Sacco, a PhD Candidate at Nanyang University of Singapore, to calculate the moduli of elasticity, shear moduli and Poisson ratios of every sample. The logic and use of the program is explained hereunder:

1. When the program is run the user can select the folder where the excel files (.csv) corresponding to all the stresses and strains are located (see Figure 8).

2. Another window is opened in which the name of headers corresponding to stress, longitudinal and axial stress must be input. By default, the names written are the headers given by the UTM for the stress and GOM Correlate for the strain respectively. It is shown in Figure 9 that ’Name of STRAIN values’ corresponds to longitudinal strain, ’Name of perpendicular STRAIN values’ to the axial strain (both obtained with GOM Correlate), and ’Name of STRESS values’ corresponds to the stress values provided by the UTM software.

(13)

(a) Select the stress folder. (b) Select the strain folder.

Figure 8: Windows to choose the folders where the excel files corresponding to stresses and strains are located.

Figure 9: The names of the STRAIN and STRESS values must coincide with the headers in the excel files.

3. Every data set of stress is coupled with its corresponding data set of strain and they are plotted one by one. In each plot the interval in which the moduli are to be calculated must be selected by chosing two stress values in the graph (Figure 10).

4. The Matlab function "robustfit" adjusts all the points selected in the previous step to a unique curve, which represents the modulus of elasticity, and its deviation. Finally, the progam records the modulus of elasticity, shear modulus and Poisson ratio of each sample and determines the average for each orientation. Therefore, running the program one time per

(14)

Figure 10: Input the lowest and highest stress value of the interval in which the moduli and Poisson ratio will be calculated.

orientation would give all the values needed to characterize the material using the equations (2-4).

Figure 11: All Stress-Strain curves for one orientation.

Figure 12: Curve of E in each sample and average E (thicker).

The value of Exis the same as E1, Ey= E2, and Ez= E3. Similarly, using equation 4 the shear moduli can be found as described in equation 4.2 of [8], and the Poisson ratios are the result of the axial or transversal elongation over the longitudinal elongation (in the direction of the strength) in the manually selected interval of the program and averaged:

ν_sample= 1 n

n

X

i=1

−ε_trans_n εlongn

(5)

(15)

where n is the number of points collected during the experiments in the selected interval.

νorientation=1 5

5

X

i=1

εn (6)

where ε_n corresponds to the Poisson ratio of each sample calculated in eq.

5.

2.2 Calculation of skin thickness

Once the material is characterized, the next step is to define the structure of the UAV. As mentioned previously, the outer design was already created, so only a new internal structure which can support the loads is studied in this section to decrease the weight of the current structure. In Table 2, some important properties and assumptions about the UAV are presented:

UAV Performance specifications

Wing span [m] 2

Payload capability [kg] 3.5

Stall angle [^◦] 15

Maximum speed [m/s] 15

Table 2: UAV performance specifications (data courtesy of Yunus Govdeli)

2.2.1 Wing geometry and determination of loads

To calculate the loads in the wing of the UAV, a model was created and ana- lyzed together with Yunus Govdeli, PhD Candidate at Nanyang University of Singapore, using the software XFLR5 with a maximum fixed speed of 15 m/s and increasing the angle of attack from 0° to 15°. XFLR5 is an analysis tool for airfoils, wings and planes operating at low Reynolds Numbers which enables a designer to estimate the aerodynamic performance of a conceptual design[9].

The geometry of the UAV is based in the already existing design in [3] but it was scaled down 50% to a wing span of 2 m (see Figure 13 and appendix 6.3).

The geometry consists on a flying wing with a Fauvel profile as shown in Figure 14, since it has a better lift performance in a wide range of Reynolds numbers;

it has a positive moment coefficient at 0° of angle of attack which improves the longitudinal stability of the plane. A comparative study with other profiles is performed in [3, p. 30].

(16)

Figure 13: Conceptual design of the UAV to be studied.

Figure 14: Side view of the drone. More detailed views of the UAV can be found in 6.3.

After simulating the UAV in XFLR5 keeping the maximum predicted speed and increasing the angle of attack from 0° to 15°, the most critical load conditions were given when the lift coefficient is highest, i.e., at 15° (stall angle of attack).

XFLR5 provides the lift and drag coefficients (Figure 15b), and bending moment along the wing span at each angle of attack every 0.5° (see Figure 15a), giving as well an approximation about drag and lift forces (Figure 15c).

(17)

(a) Bending moment along the span, in mm. (b) Lift coefficient along the span, in mm.

(c) Diagram of forces acting on the UAV, including induced drag (yellow), viscous drag (purple) and lift (green).

Figure 15: Illustrations extracted from the analysis performed with XFLR5, at 15 m/s and 15° of angle of attack.

Nonetheless, it is not possible to define in XFLR5 the internal structure of the wing. This means that the structural details had to be defined manually in different points along the span. After calculating the aerodynamic loads with XFLR5, a Matlab program was developed to estimate the net vertical forces on the structure. Loads such as the weight of electronics (battery, wires and servos), propellers and engines were input in their respective locations (one propeller in the middle of each wing and center of the plane, with the electronics). XFLR5 provides the lift coefficient for several sections along the span (76 in half of the wing in this case), which is then used to calculate the lift in each profile section following the 2-D Lift formula[10]:

L = 1

2ρC_LV_c² (7)

where ρ is the air density at the height that the aircraft is flying (1.225kg/m³), CL is the lift coefficient of a specific profile section of the span, and Vc is the relative speed of air.

(18)

Figure 16: Wing geometry and coordinate system used for the calculations.

The net vertical loads are the difference between the weight and the lift of that section along the span (axis Y in Figure 16):

Q(y) = L(y) − nW (y) (8)

W (y) is the weight and L(y) is the lift force in a particular point along the span.

n is the load factor, which can be defined as the ratio of the lift of an aircraft to its weight, and it represents conditions of the plane during a maneuver or wind gusts (n ≥ 1) or straight and level flight (n = 1). It has been selected according to [11, p. 561], which advices to be set at 2.5 for General Aviation Normal planes. It can also be calculated in accordance with its definition although the values obtained are smaller (less conservative) than 2.5.

n = L

W =

1

2ρSC_LmaxV²

W (9)

Once the forces acting on the structure have been obtained, the shear force distribution along one of the wings can be calculated integrating the net vertical force over the half span:

S(y) = Z b/2

y

q(y)dy (10)

Integrating the shear force distribution along half of the wing span, the bending moment generated is obtained:

M (y) = Z b/2

y

S(y)dy (11)

The distribution of shear force and bending moment along half of the wing span are shown in Figures 17 and 18. It is noticeable that the results given by XFLR5 (seen in Figure 15a) differ from the bending moments found with equation 11, since the weight of electronic components and inner elements, such as the spars, cannot be accurately assigned in XFLR5.

The torsional moment is calculated as descibed in [3, p. 54], ”computing the sum of the airfoil pithing moment and the moment produced due to the offset

(19)

Figure 17: Shear force distribution along half of the wing span.

Figure 18: Bending moment distribution along half of the wing.

of the elastic axis at each section from the aerodynamic center, at which the net vertical force is assumed to be acting. Thus, the moment at each cross section can be calculated as”:

Mx(y) = Cm(y)q∞c (12)

where M_x(y) is the cross-sectional pitching moment, C_mis the airfoil pitching moment coefficient given by the analysis in XFLR5, c is the chord length at that section, and q_∞ is the dynamic pressure, which represents the increase in a moving fluid’s pressure over its static value due to motion: q∞ = ^ρV₂^C²^c. Therefore, the total torsional moment, T (y) can be identified as

T (y) = M_x(y) + sQ(y) (13)

where s is the offset of the elastic axis at each section from the aerodynamic center, which is the shear center for a 2D section, assumed to be located at 25% of the chord[12]. The elastic axis is defined as a line along the span in which only bending is produced when the net force is acting on it, and it was calculated with a Matlab program designed by PhD. Candidate Yunus Govdeli [3]. The distribution is found in Figure 19:

(20)

Figure 19: Torsional moment along the span.

2.2.2 Structural sizing method

The method employed to dimension the structure shown in Figure 20 is described in [13], and it enables to obtain the skin thickness, spars thicknesses and with that information, the weight of the structure. These elements (skin, front and back spars) are modeled using one upper and one lower equivalent panels. The web of the front and back spars are also modeled separately as simple vertical panels.

Figure 20: Model of the wing box structure[13].

The equations used to determine the upper and lower skin thickness consider:

The bending moment, M, calculated previously.

Maximum allowable stresses, σmax_u and σmax_l. The maximum allowable stress in the lower panel, σmax_l, is calculated using (16) where Ex is the

(21)

modulus of elasticity and the cut off tensile strain respectively of the PLA used in orientation 3. The maximum allowable compressive stress in the upper panel, σmax_u, is calculated by multiplying the same Exand the cut off compressive strain, extracted from [13, p.107, Table 2].

The dimensionless parameter, ηtin Figure 21, relates the structural weight of the wing to the wing outer shape. The value used is 0.8, as recommended by Torenbeek[17].

The maximum thickness, tmax, of the airfoil section (see Figure 20).

The upper and lower length of the panels, Su and Sl, which can be obtained by calculating the distance between the points used to build the airfoil profile and the chord length of each section.

ts_u = M/σmax_u

ηttmaxSu

(14)

t_s_l= M/σmax_l

η_tt_maxS_l (15)

where:

σmax_l= Exεt

σmax_u= Exεc

(16)

Figure 21: Model of the wing box structure[13].

The spar thicknesses are dimensioned to withstand shear stresses. The shear flow in the walls is the sum of the shear flow due to the vertical load and the shear flow due to the torsional moment.

The flows produced by the vertical load in the front andrear spars are calculated how D. Howe proposes in Aircraft Loading and Structural Layout [18]:

q_v_{f s}= h_{f s}

h²_{f s}+ h²_rsV (17)

qv_rs = hrs

h²_{f s}+ h²_rsV (18)

where h_{f s}and h_rsare the heights of the front and rear spar respectively and V is the total vertical force obtained in (8).

The shear flow due to torsional moment is calculated as follows:

qt= T

2A_box (19)

(22)

where T is the torsional moment calculated with (13) and A_box is the area inside the wing box.

Therefore, the thickness of each spar web is calculated using the total shear flow (q = qv+ qt) and the maximum allowable shear stress, τmax:

tf s= qf s

τ_max trs= qrs

τ_max

(20)

where τ_max is the result of multiplying the laminate shear stiffness of the material in the direction that the part will be printed, G₃ and cut off shear strain, ε_s= 0.006, value extracted from [13, p.107, Table 2]:

τmax= Gxyεs (21)

At this point the thickness of the skin and spars for every particular section of the wing is known. The highest value for each web and panel will be selected to be mantained along the whole span of the wing for simplicity purposes when printing the parts. The AM process will be performed vertically (orientation 3), since according to current project members, this orientation gives the best surface finish and shortest building time. Besides, there is a minimum thickness of 0.4 mm that the printer can build, so getting a lower value would not make a big difference.

With the skin and spar thickness in the upper and lower part of each section, its width along the span and material density, given by the manufacturer, it is possible to estimate the structural weight.

W_s_u =

n

X

1

b_iS_u_it_u_maxρ_{P LA}

Ws_l=

n

X

1

biSl_itl_maxρP LA

(22)

where n is the number of sections along the span, bi is the width of the section, Su_i and Sl_i are the upper and lower panel lengths, tu_max and tl_max are the maximum upper and lower thicknesses and ρP LA is the density of PLA.

The weight of the front and rear spar is calculated as follows:

W_{f s}=

n

X

1

b_ih_{f s}_it_{f s}_maxρ_{P LA}

Wrs=

n

X

1

bihrs_itrs_maxρP LA

(23)

where tf s_max and trs_max are the maximum value of the front and rear spar thickness respectively obtained in (20), hf s_i and hrs_iare the heights of the front and rear spars.

(23)

Finally, the sum of the upper and lower spar and the skin weights is multiplied by 2 (all operations are focused on one of the wings), together with the weight of electronic components (battery, wires,etc.) and mechanical elements (propellers, control surfaces, etc.) determine the total weight of the UAV:

WT = Ws_u+ Wl_u+ Wf s+ Wrs+ We+ Wm (24)

2.3 CAD simulation and wing testing

Before building the UAV at full size, it is convenient to check the validity of the previous calculations and argumentations. For this reason, scaled wings with the same geometry and material were tested for a vertical load simulating the lift, and a 3D model has been created in SolidEdge with the real dimensions and exposed to the loads described in 2.2.1.

2.3.1 Scaled wings

Applying a vertical load on wings scaled ten times smaller than the real size (1:10) in different orientations would help the team to verify the behaviour of the new structure of the wing depending on its printing orientation and study properties related to its geometry[4].

Two samples were built in orientations 1, 2 and 3 taking as reference Figure 3, since they were the most interesting from a building point of view. Orien- tations 5 and 6 need support material, making its construction more complex and slower. Besides, they have been proved during the characterization phase to not offer more strength. Orientation 4 would not make any difference respect to orientation 1.

As the UAV is a flying wing, building this model consists on printing half of the drone stuck to a solid piece which can be held into a fixed component of the testing machine, which is the same UTM used in 2.1. The solid part of the wing is inserted in the component upside down, and another printed piece is used to push the wing downwards at 3 mm/min acting as the net force until the part breaks.

(24)

Figure 22: Scaled wing ready to be tested in the UTM.

2.3.2 A geometry issue

It appeared when testing the first batch of samples. All parts (printed only in orientation 3) broke in the sharp change of geometry in a extremely early stage of the test (see Figure 24). In a full size wing this problem could go undetected because the profile section is bigger, and a higher force would be required to produce a rupture for bending. At this small size however, only a small force is needed, given that the surface inside the section which keeps the wing stuck is smaller, and the ultimate stress, σ = F/A, is achieved with a lower force, F.

This fact is aggravated by the printing process, because when a layer is piled on another while there is a very sharp change, the filament is not properly aligned on the previous filament, making the contact surface smaller and the area which stands the pressure decrease (Figure 23). Thus, the stress increases at higher ratio as the force is increased.

Figure 23: View of filament sections piled on a straight wall on the left, and on an inclined wall at the right.

(25)

Figure 24: Scaled wings printed in orientation 3 after test.

Consequently, these geometry variations were supressed only for these tests, remaining a uniformly decreasing profile section from root to tip, and samples in orientations 1, 2 and 3 were re-built and tested with the simplified geometry (Figures 25 and 26).

Figure 25: Orientation 1 test wing with simple geometry in construction.

Figure 26: Bending of a wing with simple geometry.

2.3.3 3D Model of the UAV

A model with the real dimensions of the UAV has been designed with Solid Edge 2017 and subjected to the loads previously calculated to ensure that the structure withstand them. The model has been designed with the inner and outer structure as calculated (see Figure 27), and a new material has been added to the data base of Solid Edge with its respective properties. The simulation has been performed with tetrahedral mesh (it was necessary to simplify the model), the net load was applied in the lower face of the wing and the weight was automatically placed by the software.

(26)

Figure 27: 3D model with a section view showing the inner structure.

3 Results

The result of applying the previous theory is shown in this section. Firstly, the characterization of the material will be described and then, the result of applying that characterization into the Matlab program to obtain the thickness and weight of the UAV.

3.1 Characterization of material

Data collected during the tests was input into the Matlab program to obtain the engineering constants. However, the first time these that tests were performed, considerable differences in the stress-strain properties were seen for samples of the same orientation (see Figure 28). This could be due to a long storage time (2 weeks) in the office, exposed to high humidity (the annual average relative humidity in Singapore is 80%) and high temperature (annual average temperature is 28ºC)[19].

(27)

Figure 28: Stress-strain graphs of samples stored in the office for 2 weeks.

This difference within the results forced that the tests were repeated. This time, samples were stored inside a dry box, isolated from outer humidity, and all samples were tested within the first two days after their construction. The outcome was this time more homogeneous, as it can be seen in Figure 29.

(28)

(a) Samples printed in orientation 1. (b) Samples printed in orientation 2.

(c) Samples printed in orientation 3. (d) Samples printed in orientation 4.

(e) Samples printed in orientation 5. (f) Samples printed in orientation 6.

Figure 29: Stress-strain graphs corresponding to samples tested within the first two days after their construction, stress in MPa.

As explained in the previous section, the modulus of elasticity for each sample was deduced from a manually selected segment of the curve, normally less than 5% of the data which make up the elastic (straight) part of the slope, as it was recommended by Miquel Domingo-Espin, author of [4]. Then the average modulus for each sample was estimated, obtaining the results shown in Figure 31.

The results show that orientations 1, 2 and 4 present similar modulus of elas-

(29)

Figure 30: Average stress-strain curves for each orientation.

Figure 31: Average moduli of elasticity for each orientation.

ticity, although 4 is more fragile, as well as orientations 3 and 5, whose moduli are lower, similar to 6. Orientation 6 specimens show a similar modulus but a more plastic behaviour. The highest tensile strength is obtained with samples in orientation 2 followed by samples 4 and 1. The lowest tensile strength however is achieved when orienting the part in the 3 direction. This is related to the direction that the filaments are deposited: samples 2 show the highest value of elastic modulus and tensile strength because it has more layers pulled longi- tudinally. Samples 3 present a fragile fracture and the lowest tensile strength because the pulling direction is perpendicular to the layers, where the bonding strength is weaker than the strength of pulled contours.

It can also be deduced that the results of samples tested immediately after they have been built are more homogeneous because there are no important differences between the slopes for samples printed in the same orientation. Only some significant variations are found in the plastic behaviour of orientations 4 and 6, still for unkonwn reasons (it could be the position in which they were printed, diferences during the fixing of the piece in the testing machine, etc.), but it is irrelevant in this case because this part of the graph is not considered to calculate any modulus or Poisson ratio.

All this information was processed with the program E calculator.m for every set of samples in each orientation, and the results are collected in Table 3.

Orientation Build plane Ori.

Elastic Modulus

(MPa)

Tensile strength

(MPa)

Tensile strain (mm/mm)

Ultimate strength (MPa)

Ultimate strain (mm/mm)

Poisson ratio

1 xy x 2.85E+03 33.97 1.74E-02 28.65 1.22E-01 0.30

2 yz y 2.78E+03 38.02 1.18E-02 30.74 1.17E-01 0.25

3 xz z 2.16E+03 19.66 9.69E-03 19.62 9.73E-03 0.19

4 xy x+45 2.71E+03 33.03 1.54E-02 30.88 2.89E-02 0.26 5 yz y+45 2.28E+03 26.30 1.54E-02 24.69 1.98E-02 0.27 6 xz z+45 2.42E+03 27.52 1.63E-02 24.13 2.66E-02 0.26

Table 3: Average properties for each orientation shown in Figure 3.

To characterize the material, nine engineering constants are needed: the values of Ex, Ey and Ez are the same as E1, E2, E3 respectively. Similarly, ν1 = νxy, ν2 = νyz, and ν3 = νxz. On the other hand, the shear moduli were

(30)

obtained using (4) and orientations 4, 5 and 6 for G_xy, G_yzand G_xzrespectively.

Ex 2.85 · 10³ νxy 0.30 Gxy 1.08 · 10³ Ey 2.78 · 10³ νyz 0.25 Gyz 9.61 · 10² Ez 2.16 · 10³ νxz 0.19 Gxz 8.96 · 10²

Table 4: Engineering constants

These results represent the shear and elastic moduli, and Poisson ratios in each direction according to Figure 3. Nevertheless, as it has been mentioned in the methodology section (2.2.1), calculations were based on z-direction and the xz - plane (orientation 3), since most of the pieces will be built in that orientation for complexity and functional reasons also described before. Besides, this choice is conservative, since even if there was a part built in a different orientation, it would present stronger mechanical properties.

3.2 Structural sizing

Following the steps described in 2.2.2, the final skin and spar thickness will be the maximum value along the span (see Figures 32 and 33). It can be observed that maximum skin thickness is achieved at 0.1 m away from the axis of symmetry of the UAV, where the sharp change of geometry occurs. This thickness is slightly higher in the upper face due to the geometry of the wing;

the upper panel is shorter than the lower panel. This maximum value is 0.51 mm, which will be approximated to 0.60 mm so it can be printed, and it will be the skin thickness in the whole wing and both faces.

The maximum spar thickness is given at the front spar, located where the airfoil is thickest (about 25% of the chord in this case). The maximum value was given near the ”fuselage”, where the electronic components will be placed and where the abrupt geometry change is located. This value was 0.81 mm, which was rounded to 1.00 mm, a more conservative and feasible thickness to build with a printer. On the other hand, the rear spar was located at about two thirds from the leading edge, and this thickness resulted equal to 0.53 mm, which was rounded to 0.60 mm.

Figure 32: Variation of the skin thickness along half of the span.

Figure 33: Variation of the spar thickness along half of the span.

(31)

The values were calculated with the code shown in Appendix 6.2, and are summarized in the following table:

Obtained Rectified

Skin thickness (mm) 0.51 0.60

Spar thickness(mm) Front spar 0.81 1.00

Rear spar 0.53 0.60

Structural weight (kg) 0.89 1.05

Total weight (kg) 2.97 3.13

Table 5: Obtained values by the Matlab program and modified values to facilitate the printing process.

Values classified as Rectified in Table 5 are the obtained values rounded to the next higher closest value (more conservative) which can be produced by the printer.

The structural weight calculated with the rectified values is 0.5 kg lighter than the prototype already existing [3].

3.3 Simulation

The first set of three scaled wings built in orientation 3 (”Orientation 3 Com- plex”) presented some geometrical issues due to abrupt variations in the profile near to the axis of symmetry of the UAV not only during the test but also during the printing process. The geometry was then simplified, built (Figure 34) and tested in orientations 1, 2 and 3. Two samples of each orientation were tested and their results averaged. Values are compared in Figure 35 and Table 6.

(a) Orientation 1. (b) Orientation 2. (c) Orientation 3.

Figure 34: Orientation of the printed parts.

In Figure 35, the average values obtained during the tests are collected. Ori- entation 3 C corresponds to the first samples tested with complex geometry, and it can be seen that there is a noticeable difference in the wing strength between them and Orientation 3 S, the simplified samples, although the maximum displacement that they admit is similar. In the simplified structures, Orientation 1 presents the highest strength and diplacement, but similar stiffness to Orienta- tion 3 before they break. Orientation 1 is also the most elastic and Orientation 3 the most plastic.

(32)

Figure 35: Average behaviour in bending test of scaled wings. The letter ”S”

stands for the wing with simplified geometry whilst ”C” for the original shape.

The orientation which will probably be used the most during the construction of the UAV will be Orientation 3 for its simplicity and building speed, and according to Figure 35 and Table 6, this orientation would stand forces almost as high as orientations 1 and 2, although the break would arrive at smaller displacements while the wing is bending.

Sample Fmax_avg (N) dmax_avg (mm)

Orientation 1 46.02 37.42

Orientation 3 (Complex) 5.44 9.24

Table 6: Averaged maximum force and displacement of each orientation

3.3.1 3D Model

The purpose of this simulation is to verify that the stress on the surface does not achieve the ultimate strength and neither the yield strength (after this value the material deforms plastically) given by the properties determined during the characterization. The simulation considers the vertical net forces (Figure 36), which determine the skin thickness.

(33)

Figure 36: Net forces considered for the simulation.

Lift forces (84 N) are distributed at the bottom face of the wings while the weight (31.3 N) is automatically placed in the center of gravity and the fuselage, where the electronic components are placed, has been fixed.

The result can be seen in Figure 37; the maximum stress is achieved at the wings roots, and it decreases spanwise towards the tips. However, this maximum stress does not exceed the ultimate and yield stress. The inner geometry has been simplified, deleting the spars to enable the meshing and hence, the results may vary from reality, but it is still appreciable that the maximum stress in the surface is considerably lower than the yield stress (19.6 MPa).

Figure 37: Stress on the UAV surface simulating a flight at 15 m/s and 15° of angle of attack in Solid Edge.

(34)

4 Conclusion

The work performed in this report has been quite theoretical excepting the characterization part. However, some details can be inferred from this dissertation.

Characterization In the first phase, it was deduced that the material is vir- tually isotropic although deviations in the results depending on the orientation appeared due to uneven filling and deposition during the printing process.

During testing, samples with the thinnest gauge (D638 - Type IV) were proved to provide the same results as bigger samples, which take a longer construction time, and the rupture is easier to detect by camera.

The impact of external conditions on the material must also be stressed; samples were first tested after they had been stored in an office for two or three weeks, and the results showed considerable irregularities compared with the same samples tested within a maximum of two days after manufacturing. This means that the structure of the UAV could be compromised after long exposure to high temperature and/or moist, and it is advisable to either look for a material less sensitive to these external factors or to more throughly investigate these ef- fects, since they are likely to strongly affect both the stiffness and the strength of the entire aircraft.

Structural sizing and simulation Taking the properties of the theoretically lowest strength orientation, the weight of the UAV has been decreased with this sizing method more than half a kilogram compared with previous structures studied in the project[3]. It is till possible to perform more complete analysis and a more faithful representation of the model but these calculations are promising.

On the other hand, the tests of the scaled wings revealed that a wing would be significantly stronger if the abrupt changes of geometry were suppressed and the wing root decreased uniformly from the root to the tip. This fact was aggravated in the scaled model since the size was much smaller and this could possibly influence in the printing process. This is not an easy task, because the box in which all the electronic components are, does not have a complete profile (it would mean more weight in that case), but this could lead to start to study a smoother transition along the span.

(35)

5 Future work

Time is always a limiting a factor, and it has not been an exception in this work.

Thus, these are some aspects which can be developed further on:

1. Automate the Matlab program E calculator.m to a higher extent: right now, the section in which the moduli and Poisson ratios are calculated is selected manually. One idea to develop in this respect would be to im- plement an algorythm to detect slope variations between following points and make the selected section stop in that point.

2. Influence of humidity and/or temperature on the UAV structure: according to the results during characterization, storing the samples in an office gave irregular results compared with the ones stored in a dry box and tested within two days. Therefore, it is recommended to further study how these factors are affecting the material and the structure, or consider another material more resistant to these conditions.

3. A more complete simulation.

(a) More loads can be taken into consideration, such as drag.

(b) The location of loads can be more precise once the exact weight and position of components (battery, wires, propellers, etc.) are known.

4. Construction of the UAV: all the work executed and described in this report has had the final purpose of building a full size functional structure.

Therefore, this point is to be performed.

(36)

6 Appendix

6.1 E calculator.m

1 % Startup

2 clear all;

3 close all;

4

5 disp('Select STRESS folder');

6 stress folder = uigetdir(path,'Select STRESS folder');

7 disp('Select STRAIN folder');

8 strain folder = uigetdir(path,'Select STRAIN folder');

9

10 prompt = {'Plot name:','x label (strain):','y label ...

(stress):','Name of STRAIN values:','Name of STRESS ...

values:','Name of perpendicular STRAIN values:'};

11 dlg title = 'Input';

12 num lines = 1;

13 defaultans = {'Stress vs Strain','Strain','Stress','Surface ...

component 1.avg(epsY) [True strain]','Stress','Surface ...

component 1.avg(epsX) [True strain]'};

14 answer = inputdlg(prompt,dlg title,num lines,defaultans);

15

16 filePattern = sprintf('%s/*.csv', stress folder);

17 stress file list = dir(filePattern);

18 filePattern = sprintf('%s/*.csv', strain folder);

19 strain file list = dir(filePattern);

20

21 name plot = answer{1};

22 x label = answer{2};

23 y label = answer{3};

24 strain col = answer{4};

25 stress col = answer{5};

26 strainX col = answer{6};

27

28 des cols = {stress col,strain col};

29 smallest n = 100000;

30 biggest n = 1;

31 avg stress values = [];

32 avg strain values = [];

33

34 % Determine the biggest and smallest set of data

35 for file = stress file list'

36 M = xlsread(fullfile(stress folder,file.name));

37 [row,col] = size(M);

38 if smallest n > row

39 smallest n = row;

40 end

41 if biggest n< row

42 biggest n = row;

43 end

44 end

45 for file = strain file list'

46 M=xlsread(fullfile(strain folder,file.name));

47 [row,col]=size(M);

48 if smallest n > row

49 smallest n = row;

50 end

51 if biggest n< row

52 biggest n = row;

(37)

53 end

54 end

55

56 stress data = cell(biggest n,length(strain file list));

57 strain data = cell(biggest n,length(strain file list));

58 strainX data = cell(biggest n,length(strain file list));

59

60 % Extract columns

61 ii=0;

62 for file = stress file list'

63 ii = ii + 1;

64 [n,s,r] = xlsread(fullfile(stress folder,file.name));

65 [row, col] = find(strcmpi(s,stress col));

66 stress values = n(1:end,col);

67 avg stress values(:,ii)=stress values(1:smallest n);

68 columns to pad = biggest n - length(stress values);

69 padding = num2cell(NaN*ones(columns to pad,1));

70 stress values = num2cell(stress values);

71 padded col = [stress values; padding];

72 stress data(:,ii) = padded col;

73 end

74 ii=0;

76 ii = ii + 1;

77 [n,s,r] = xlsread(fullfile(strain folder,file.name));

78 [row, col] = find(strcmpi(s,strain col));

79 strain values = n(1:end,col);

80 avg strain values(:,ii)=strain values(1:smallest n);

81 columns to pad = biggest n - length(strain values);

83 strain values = num2cell(strain values);

84 padded col = [strain values; padding];

85 strain data(:,ii) = padded col;

86

87 [row, col] = find(strcmpi(s,strainX col));

88 strainX values = n(1:end,col);

89 avg strainX values(:,ii)=strainX values(1:smallest n);

91 strainX values = num2cell(strainX values);

92 padded col = [strainX values; padding];

93 strainX data(:,ii) = padded col;

94 end

95

96 % Plot data of each test

97 98 ii=0;

99 strain data = cell2mat(strain data);

100 strainX data = cell2mat(strainX data);

101 stress data = cell2mat(stress data);

102 new stress data = [];

103 new strain data = [];

104 sum ny = [];

105 ny = [];

107 ii = ii + 1;

108 plot(strain data(:,ii),stress data(:,ii),'DisplayName',file.name);

109 grid on;

110 legend('Location','SE')

111 set(gca,'FontSize',20);

112 ylabel({y label});

113 xlabel({x label});

114 title({name plot});

(38)

115

116 % pause;

117

118 % Point elimination

119 prompt = {'Eliminate values above this STRESS:','Eliminate ...

values below this STRESS:',};

120 dlg title = 'Do you want to eliminate values bigger and ...

smaller than a certain stress?';

121 num lines = 1;

122 defaultans = {'0','0'};

123 answer2 = inputdlg(prompt,dlg title,num lines,defaultans);

124

125 stress points to keepB = 1000000;

126 stress points to keepS = 1;

127 strain points to elim = 1000000;

128

129 col = 0;

130 row = 0;

131 [rows,cols]=size(stress data);

132 % Eliminate points using STRESS

133 if answer2{1} 6= 0

134 col = ii;

135 ixB = find(stress data(:,ii)> ...

str2num(answer2{1}),1,'first');

136 if stress points to keepB > ixB

137 stress points to keepB = ixB;

138 end

139 end

140

141 if answer2{2} 6= 0

142 col = ii;

143 ixS = find(stress data(:,ii)> ...

str2num(answer2{2}),1,'first');

144 if ixS ≤ min(stress data(:,ii))

145 ixS = 1;

146 end

147 end

148 sumStress = stress data(ixS:ixB,ii);

149 sumStrain = strain data(ixS:ixB,ii);

150 sumStrainX = strainX data(ixS:ixB,ii);

151 sum ny = -(sumStrainX./sumStrain);

152 ny(1,ii) = nanmean(sum ny);

153 columns to pad = length(stress data) - length(sumStress);

154 padding = NaN*ones(columns to pad,1);

155 new stress data(:,ii) = [sumStress; padding];

156 new strain data(:,ii) = [sumStrain; padding];

157 new ny(:,ii) = [sum ny;padding];

158 159 end

160 avg ny = mean(ny);

161 close all;

162 ii = 0;

163 j=0;

164 ultStrain = []; ultStress = []; TStrain = []; TStress = [];

165 f1 = figure;

166 C = {'k','b','r','g','y',[.5 .6 .7],[.8 .2 .6]}; % Cell array of colros

167 legendtext = {'Sample 1','Sample 2','Sample 3','Sample 4','Sample ...

5','Sample 6','Sample 7','Sample 8'};

169 ii = ii + 1;

170 h1(ii) = plot(strain data(:,ii),stress data(:,ii),'color',C{ii});

171 grid on;

(39)

172 legend(h1,legendtext(1:ii))

173 hold on;

178

179 [ultStrain(1,ii),j] = max(strain data(:,ii));

180 ultStress(1,ii) = (stress data(j,ii));

181

182 [TStress(1,ii),j] = max(stress data(:,ii));

183 TStrain(1,ii) = strain data(j,ii);

184 185 end

186

187 f2 = figure;

188 Legend = cell(ii,1);

189 ii = 0;

190 E = [];

191 y0 = [];

192 [maxim,i] = max(max(new strain data));

193 strain plot = (0:maxim/50:maxim*1.5);

194 195 196

198 ii = ii + 1;

199 scatter(new strain data(:,ii),new stress data(:,ii),[],C{ii},'filled');

200 [brob,stats] = ...

robustfit(new strain data(:,ii),new stress data(:,ii));

201 grid on; hold on;

202 h2(ii) = ...

plot(strain plot,brob(1)+brob(2)*strain plot,'color',C{ii},'LineWidth',1);

203 legend(h2,legendtext(1:ii))

204 y0(ii) = brob(1);

205 E(ii) = brob(2);

206 end

207

208 avg TStress = mean(TStress);

209 avg TStrain = mean(TStrain);

210 avg ultStress = mean(ultStress);

211 avg ultStrain = mean(ultStrain);

212 avg E = mean(E);

213 avg y0 = mean(y0);

214 plot(strain plot,avg y0 + avg E * strain plot,'LineWidth',4);

215 hold on;

6.2 Thickness calculation

6.2.1 Main.m

1 clear all;

2 %% Reading profile of wing

3 Fauvel = xlsread('Fauvel.xlsx');

4 FauX = Fauvel (1:end,1); % The X coordinates of the Fauvel ...

wing are not the same in the upper and lower part

5 FauY = Fauvel (1:end,2);

Characterization of PLA and design of a 3D printed wing