First commit!

Changed coordinate lists into numpy arrays.

Commit
588c34a3d595fcad5e93b8d4893f1098ce64d046
Author
blendoit <blendoit@gmail.com>
Author date
Committer
blendoit <blendoit@gmail.com>
Committer date
Changed files
.gitignore
index 00000000..c7969028 000000..100644
@@ -0,0 +1,4 @@
1 Added: # .gitignore
2 Added: **/__pycache__/
3 Added: **/log.txt
4 Added: save/
README.org
index 00000000..caa04a98 000000..100644
@@ -0,0 +1,7 @@
1 Added: #+TITLE: UCLA MAE 154B
2 Added: #+SUBTITLE: Spring 2019 Final Project
3 Added:
4 Added: This program enables the creation of NACA airfoils;
5 Added: the analysis of the airfoil's structural properties;
6 Added: the optimization via genetic algorithm of a population of airfoils;
7 Added: With the final objective of designing a lightweight FAR 23 compliant airfoil.
creator/fuselage.py
index 00000000..e69de29b 000000..100644
creator/propulsion.py
index 00000000..e69de29b 000000..100644
creator/wing.py
index 00000000..4988cb5d 000000..100644
@@ -0,0 +1,375 @@
1 Added: """
2 Added: The wing.py module contains class definitions for and various components
3 Added: we add to an airfoil (spars, stringers, and ribs).
4 Added:
5 Added: Classes:
6 Added: Airfoil: instantiated with class method to provide coordinates to heirs.
7 Added: Spar: inherits from Airfoil.
8 Added: Stringer: also inherits from Airfoil.
9 Added:
10 Added: Functions:
11 Added: plot_geom(airfoil): generates a 2D plot of the airfoil & any components.
12 Added: """
13 Added:
14 Added: import sys
15 Added: import os.path
16 Added: import logging
17 Added: import numpy as np
18 Added: from math import sin, cos, atan
19 Added: import bisect as bi
20 Added: import matplotlib.pyplot as plt
21 Added:
22 Added: logging.basicConfig(filename='log.txt',
23 Added: level=logging.DEBUG,
24 Added: format='%(asctime)s - %(levelname)s - %(message)s')
25 Added:
26 Added:
27 Added: class Component:
28 Added: """Basic component providing coordinates and tools."""
29 Added:
30 Added: # TODO: define defaults in separate module
31 Added: def __init__(self):
32 Added: self.x = np.array([])
33 Added: self.z = np.array([])
34 Added: self.material = str()
35 Added: self.mass = float()
36 Added:
37 Added: def set_material(self, material):
38 Added: """Set the component bulk material."""
39 Added: self.material = material
40 Added:
41 Added: def info_print(self, round):
42 Added: """Print all the component's coordinates to the terminal."""
43 Added: name = f' CREATOR DATA FOR {str(self).upper()} '
44 Added: num_of_dashes = len(name)
45 Added: print(num_of_dashes * '-')
46 Added: print(name)
47 Added: for k, v in self.__dict__.items():
48 Added: if type(v) != list:
49 Added: print('{}:\n'.format(k), v)
50 Added: print(num_of_dashes * '-')
51 Added: for k, v in self.__dict__.items():
52 Added: if type(v) == list:
53 Added: print('{}:\n'.format(k), np.around(v, round))
54 Added: return None
55 Added:
56 Added: def info_save(self, save_path, number):
57 Added: """Save all the object's coordinates (must be full path)."""
58 Added: file_name = f'{str(self).lower()}_{number}.txt'
59 Added: full_path = os.path.join(save_path, file_name)
60 Added: try:
61 Added: with open(full_path, 'w') as sys.stdout:
62 Added: self.info_print(6)
63 Added: # This line required to reset behavior of sys.stdout
64 Added: sys.stdout = sys.__stdout__
65 Added: logging.debug(f'Successfully wrote to file {full_path}')
66 Added: except IOError:
67 Added: print(f'Unable to write {file_name} to specified directory.\n',
68 Added: 'Was the full path passed to the function?')
69 Added: return None
70 Added:
71 Added:
72 Added: class Airfoil(Component):
73 Added: """This class represents a single NACA airfoil.
74 Added:
75 Added: The coordinates are saved as two lists
76 Added: for the x- and z-coordinates. The coordinates start at
77 Added: the leading edge, travel over the airfoil's upper edge,
78 Added: then loop back to the leading edge via the lower edge.
79 Added:
80 Added: This method was chosen for easier future exports
81 Added: to 3D CAD packages like SolidWorks, which can import such
82 Added: geometry as coordinates written in a CSV file.
83 Added: """
84 Added:
85 Added: # TODO: default values in separate module
86 Added: def __init__(self, chord, semi_span, material):
87 Added: super().__init__()
88 Added: # self.x = np.array([])
89 Added: # self.z = np.array([])
90 Added: # self.chord = chord
91 Added: """Create airfoil from its chord and semi-span."""
92 Added: self.chord = chord if chord > 20 else 20
93 Added: if chord <= 20:
94 Added: logging.debug('Chord too small, using minimum value of 20.')
95 Added: self.semi_span = semi_span
96 Added: self.material = material
97 Added: self.naca_num = int()
98 Added:
99 Added: def __str__(self):
100 Added: return type(self).__name__
101 Added:
102 Added: def add_naca(self, naca_num):
103 Added: """Generate surface geometry for a NACA airfoil.
104 Added:
105 Added: The nested functions perform the required steps to generate geometry,
106 Added: and can be called to solve the geometry y-coordinate for any 'x' input.
107 Added: Equation coefficients were retrieved from Wikipedia.org.
108 Added:
109 Added: Parameters:
110 Added: naca_num: 4-digit NACA wing
111 Added:
112 Added: Return:
113 Added: None
114 Added: """
115 Added: self.naca_num = naca_num
116 Added: # Variables extracted from naca_num argument passed to the function
117 Added: m = int(str(naca_num)[0]) / 100
118 Added: p = int(str(naca_num)[1]) / 10
119 Added: t = int(str(naca_num)[2:]) / 100
120 Added: # x-coordinate of maximum camber
121 Added: p_c = p * self.chord
122 Added:
123 Added: def get_camber(x):
124 Added: """
125 Added: Returns camber z-coordinate from 1 'x' along the airfoil chord.
126 Added: """
127 Added: z_c = float()
128 Added: if 0 <= x < p_c:
129 Added: z_c = (m / (p**2)) * (2 * p * (x / self.chord) -
130 Added: (x / self.chord)**2)
131 Added: elif p_c <= x <= self.chord:
132 Added: z_c = (m /
133 Added: ((1 - p)**2)) * ((1 - 2 * p) + 2 * p *
134 Added: (x / self.chord) - (x / self.chord)**2)
135 Added: return (z_c * self.chord)
136 Added:
137 Added: def get_thickness(x):
138 Added: """Return thickness from 1 'x' along the airfoil chord."""
139 Added: x = 0 if x < 0 else x
140 Added: z_t = 5 * t * self.chord * (+0.2969 *
141 Added: (x / self.chord)**0.5 - 0.1260 *
142 Added: (x / self.chord)**1 - 0.3516 *
143 Added: (x / self.chord)**2 + 0.2843 *
144 Added: (x / self.chord)**3 - 0.1015 *
145 Added: (x / self.chord)**4)
146 Added: return z_t
147 Added:
148 Added: def get_theta(x):
149 Added: dz_c = float()
150 Added: if 0 <= x < p_c:
151 Added: dz_c = ((2 * m) / p**2) * (p - x / self.chord)
152 Added: elif p_c <= x <= self.chord:
153 Added: dz_c = (2 * m) / ((1 - p)**2) * (p - x / self.chord)
154 Added:
155 Added: theta = atan(dz_c)
156 Added: return theta
157 Added:
158 Added: def get_coord_u(x):
159 Added: x = x - get_thickness(x) * sin(get_theta(x))
160 Added: z = get_camber(x) + get_thickness(x) * cos(get_theta(x))
161 Added: return (x, z)
162 Added:
163 Added: def get_coord_l(x):
164 Added: x = x + get_thickness(x) * sin(get_theta(x))
165 Added: z = get_camber(x) - get_thickness(x) * cos(get_theta(x))
166 Added: return (x, z)
167 Added:
168 Added: # Densify x-coordinates 10 times for first 1/4 chord length
169 Added: x_chord_25_percent = round(self.chord / 4)
170 Added: x_chord = [i / 10 for i in range(x_chord_25_percent * 10)]
171 Added: x_chord.extend(i for i in range(x_chord_25_percent, self.chord + 1))
172 Added: # Generate our airfoil skin geometry from previous sub-functions
173 Added: self.x_c = np.array([])
174 Added: self.z_c = np.array([])
175 Added: # Upper surface and camber line
176 Added: for x in x_chord:
177 Added: self.x_c = np.append(self.x_c, x)
178 Added: self.z_c = np.append(self.z_c, get_camber(x))
179 Added: self.x = np.append(self.x, get_coord_u(x)[0])
180 Added: self.z = np.append(self.z, get_coord_u(x)[1])
181 Added: # Lower surface
182 Added: for x in x_chord[::-1]:
183 Added: self.x = np.append(self.x, get_coord_l(x)[0])
184 Added: self.z = np.append(self.z, get_coord_l(x)[1])
185 Added: return None
186 Added:
187 Added:
188 Added: class Spar(Component):
189 Added: """Contains a single spar's data."""
190 Added: def __init__(self, airfoil, loc_percent, material):
191 Added: """Set spar location as percent of chord length."""
192 Added: super().__init__()
193 Added: super().set_material(material)
194 Added: self.cap_area = float()
195 Added: loc = loc_percent * airfoil.chord
196 Added: # bi.bisect_left: returns index of first value in airfoil.x > loc
197 Added: # This ensures that spar geom intersects with airfoil geom.
198 Added: # Spar upper coordinates
199 Added: spar_u = bi.bisect_left(airfoil.x, loc) - 1
200 Added: self.x = np.append(self.x, airfoil.x[spar_u])
201 Added: self.z = np.append(self.z, airfoil.z[spar_u])
202 Added: # Spar lower coordinates
203 Added: spar_l = bi.bisect_left(airfoil.x[::-1], loc)
204 Added: self.x = np.append(self.x, airfoil.x[-spar_l])
205 Added: self.z = np.append(self.z, airfoil.z[-spar_l])
206 Added: return None
207 Added:
208 Added: def set_cap_area(self, cap_area):
209 Added: self.cap_area = cap_area
210 Added: return None
211 Added:
212 Added: def set_mass(self, mass):
213 Added: self.mass = mass
214 Added: return None
215 Added:
216 Added:
217 Added: class Stringer(Component):
218 Added: """Contains the coordinates of all stringers."""
219 Added: def __init__(self):
220 Added: super().__init__()
221 Added: self.x_start = []
222 Added: self.x_end = []
223 Added: self.z_start = []
224 Added: self.z_end = []
225 Added: self.diameter = float()
226 Added: self.area = float()
227 Added:
228 Added: def add_coord(self, airfoil, spars, stringer_u_1, stringer_u_2,
229 Added: stringer_l_1, stringer_l_2):
230 Added: """Add equally distributed stringers to four airfoil locations
231 Added: (upper nose, lower nose, upper surface, lower surface).
232 Added:
233 Added: Parameters:
234 Added: airfoil_coord: packed airfoil coordinates
235 Added: spar_coord: packed spar coordinates
236 Added: stringer_u_1: upper nose number of stringers
237 Added: stringer_u_2: upper surface number of stringers
238 Added: stringer_l_1: lower nose number of stringers
239 Added: stringer_l_2: lower surface number of stringers
240 Added:
241 Added: Returns:
242 Added: None
243 Added: """
244 Added:
245 Added: # Find distance between leading edge and first upper stringer
246 Added: interval = spars.x[0][0] / (stringer_u_1 + 1)
247 Added: # initialise first self.stringer_x at first interval
248 Added: x = interval
249 Added: # Add upper stringers from leading edge until first spar.
250 Added: for _ in range(0, stringer_u_1):
251 Added: # Index of the first value of airfoil.x > x
252 Added: i = bi.bisect_left(airfoil.x, x)
253 Added: self.x.append(airfoil.x[i])
254 Added: self.z.append(airfoil.z[i])
255 Added: x += interval
256 Added: # Add upper stringers from first spar until last spar
257 Added: # TODO: stringer placement if only one spar is created
258 Added: interval = (airfoil.spar.x[-1][0] -
259 Added: airfoil.spar.x[0][0]) / (stringer_u_2 + 1)
260 Added: x = interval + airfoil.spar.x[0][0]
261 Added: for _ in range(0, stringer_u_2):
262 Added: i = bi.bisect_left(airfoil.x, x)
263 Added: self.x.append(airfoil.x[i])
264 Added: self.z.append(airfoil.z[i])
265 Added: x += interval
266 Added:
267 Added: # Find distance between leading edge and first lower stringer
268 Added: interval = airfoil.spar.x[0][1] / (stringer_l_1 + 1)
269 Added: x = interval
270 Added: # Add lower stringers from leading edge until first spar.
271 Added: for _ in range(0, stringer_l_1):
272 Added: i = bi.bisect_left(airfoil.x[::-1], x)
273 Added: self.x.append(airfoil.x[-i])
274 Added: self.z.append(airfoil.z[-i])
275 Added: x += interval
276 Added: # Add lower stringers from first spar until last spar
277 Added: interval = (airfoil.spar.x[-1][1] -
278 Added: airfoil.spar.x[0][1]) / (stringer_l_2 + 1)
279 Added: x = interval + airfoil.spar.x[0][1]
280 Added: for _ in range(0, stringer_l_2):
281 Added: i = bi.bisect_left(airfoil.x[::-1], x)
282 Added: self.x.append(airfoil.x[-i])
283 Added: self.z.append(airfoil.z[-i])
284 Added: x += interval
285 Added: return None
286 Added:
287 Added: def add_area(self, area):
288 Added: self.area = area
289 Added: return None
290 Added:
291 Added: def add_mass(self, mass):
292 Added: self.mass = len(self.x) * mass + len(self.x) * mass
293 Added: return None
294 Added:
295 Added: def add_webs(self, thickness):
296 Added: """Add webs to stringers."""
297 Added: for _ in range(len(self.x) // 2):
298 Added: self.x_start.append(self.x[_])
299 Added: self.x_end.append(self.x[_ + 1])
300 Added: self.z_start.append(self.z[_])
301 Added: self.z_end.append(self.z[_ + 1])
302 Added: self.thickness = thickness
303 Added: return None
304 Added:
305 Added: def info_print(self, round):
306 Added: super().info_print(round)
307 Added: print('Stringer Area:\n', np.around(self.area, round))
308 Added: return None
309 Added:
310 Added:
311 Added: def plot_geom(airfoil, spars, stringers):
312 Added: """This function plots the airfoil's + sub-components' geometry."""
313 Added: fig, ax = plt.subplots()
314 Added:
315 Added: # Plot chord
316 Added: x = [0, airfoil.chord]
317 Added: y = [0, 0]
318 Added: ax.plot(x, y, linewidth='1')
319 Added: # Plot quarter chord
320 Added: ax.plot(airfoil.chord / 4,
321 Added: 0,
322 Added: '.',
323 Added: color='g',
324 Added: markersize=24,
325 Added: label='Quarter-chord')
326 Added: # Plot mean camber line
327 Added: ax.plot(airfoil.x_c,
328 Added: airfoil.z_c,
329 Added: '-.',
330 Added: color='r',
331 Added: linewidth='2',
332 Added: label='Mean camber line')
333 Added: # Plot airfoil surfaces
334 Added: ax.plot(airfoil.x, airfoil.z, color='b', linewidth='1')
335 Added:
336 Added: # Plot spars
337 Added: try:
338 Added: for spar in spars:
339 Added: x = (spar.x)
340 Added: y = (spar.z)
341 Added: ax.plot(x, y, '-', color='y', linewidth='4')
342 Added: except AttributeError:
343 Added: print('No spars to plot.')
344 Added: # Plot stringers
345 Added: try:
346 Added: for _ in range(0, len(airfoil.stringer.x)):
347 Added: x = airfoil.stringer.x[_]
348 Added: y = airfoil.stringer.z[_]
349 Added: ax.plot(x, y, '.', color='y', markersize=12)
350 Added: except AttributeError:
351 Added: print('No stringers to plot.')
352 Added:
353 Added: # Graph formatting
354 Added: # plot_bound = np.amax(airfoil.x)
355 Added: ax.set(
356 Added: title='NACA ' + str(airfoil.naca_num) + ' airfoil',
357 Added: xlabel='X axis',
358 Added: # xlim=[-0.10 * plot_bound, 1.10 * plot_bound],
359 Added: ylabel='Z axis')
360 Added: # ylim=[-(1.10 * plot_bound / 2), (1.10 * plot_bound / 2)])
361 Added:
362 Added: plt.grid(axis='both', linestyle=':', linewidth=1)
363 Added: plt.gca().set_aspect('equal', adjustable='box')
364 Added: plt.gca().legend(bbox_to_anchor=(1, 1),
365 Added: bbox_transform=plt.gcf().transFigure)
366 Added: plt.show()
367 Added: return fig, ax
368 Added:
369 Added:
370 Added: def main():
371 Added: return None
372 Added:
373 Added:
374 Added: if __name__ == '__main__':
375 Added: main()
evaluator.py
index 00000000..85325a95 000000..100644
@@ -0,0 +1,284 @@
1 Added: """
2 Added: The evaluator.py module contains a single Evaluator class,
3 Added: which knows all the attributes of a specified Airfoil instance,
4 Added: and contains functions to analyse the airfoil's geometrical
5 Added: & structural properties.
6 Added: """
7 Added:
8 Added: import sys
9 Added: import os.path
10 Added: import numpy as np
11 Added: from math import sqrt
12 Added: import matplotlib.pyplot as plt
13 Added:
14 Added:
15 Added: class Evaluator:
16 Added: """Performs structural evaluations for the airfoil passed as argument."""
17 Added: def __init__(self, airfoil):
18 Added: # Evaluator knows all geometrical info from evaluated airfoil
19 Added: self.airfoil = airfoil
20 Added: self.spar = airfoil.spar
21 Added: self.stringer = airfoil.stringer
22 Added: # Global dimensions
23 Added: self.chord = airfoil.chord
24 Added: self.semi_span = airfoil.semi_span
25 Added: # Mass & spanwise distribution
26 Added: self.mass_total = float(airfoil.mass + airfoil.spar.mass +
27 Added: airfoil.stringer.mass)
28 Added: self.mass_dist = []
29 Added: # Lift
30 Added: self.lift_rectangular = []
31 Added: self.lift_elliptical = []
32 Added: self.lift_total = []
33 Added: # Drag
34 Added: self.drag = []
35 Added: # centroid
36 Added: self.centroid = []
37 Added: # Inertia terms:
38 Added: self.I_ = {'x': 0, 'z': 0, 'xz': 0}
39 Added:
40 Added: def __str__(self):
41 Added: return type(self).__name__
42 Added:
43 Added: def info_print(self, round):
44 Added: """Print all the component's evaluated data to the terminal."""
45 Added: name = ' EVALUATOR DATA FOR {} '.format(str(self).upper())
46 Added: num_of_dashes = len(name)
47 Added: print(num_of_dashes * '-')
48 Added: print(name)
49 Added: for k, v in self.__dict__.items():
50 Added: if type(v) != list:
51 Added: print('{}:\n'.format(k), v)
52 Added: print(num_of_dashes * '-')
53 Added: for k, v in self.__dict__.items():
54 Added: if type(v) == list:
55 Added: print('{}:\n'.format(k), np.around(v, round))
56 Added: return None
57 Added:
58 Added: def info_save(self, save_path, number):
59 Added: """Save all the object's coordinates (must be full path)."""
60 Added: file_name = 'airfoil_{}_eval.txt'.format(number)
61 Added: full_path = os.path.join(save_path, file_name)
62 Added: try:
63 Added: with open(full_path, 'w') as sys.stdout:
64 Added: self.info_print(6)
65 Added: # This line required to reset behavior of sys.stdout
66 Added: sys.stdout = sys.__stdout__
67 Added: print('Successfully wrote to file {}'.format(full_path))
68 Added: except IOError:
69 Added: print(
70 Added: 'Unable to write {} to specified directory.\n'.format(
71 Added: file_name), 'Was the full path passed to the function?')
72 Added: return None
73 Added:
74 Added: # All these functions take integer arguments and return lists.
75 Added:
76 Added: def get_lift_rectangular(self, lift):
77 Added: L_prime = [lift / (self.semi_span * 2) for x in range(self.semi_span)]
78 Added: return L_prime
79 Added:
80 Added: def get_lift_elliptical(self, L_0):
81 Added: L_prime = [
82 Added: L_0 / (self.semi_span * 2) * sqrt(1 - (y / self.semi_span)**2)
83 Added: for y in range(self.semi_span)
84 Added: ]
85 Added: return L_prime
86 Added:
87 Added: def get_lift_total(self):
88 Added: F_z = [(self.lift_rectangular[_] + self.lift_elliptical[_]) / 2
89 Added: for _ in range(len(self.lift_rectangular))]
90 Added: return F_z
91 Added:
92 Added: def get_mass_distribution(self, total_mass):
93 Added: F_z = [total_mass / self.semi_span for x in range(0, self.semi_span)]
94 Added: return F_z
95 Added:
96 Added: def get_drag(self, drag):
97 Added: # Transform semi-span integer into list
98 Added: semi_span = [x for x in range(0, self.semi_span)]
99 Added:
100 Added: # Drag increases after 80% of the semi_span
101 Added: cutoff = round(0.8 * self.semi_span)
102 Added:
103 Added: # Drag increases by 25% after 80% of the semi_span
104 Added: F_x = [drag for x in semi_span[0:cutoff]]
105 Added: F_x.extend([1.25 * drag for x in semi_span[cutoff:]])
106 Added: return F_x
107 Added:
108 Added: def get_centroid(self):
109 Added: """Return the coordinates of the centroid."""
110 Added: stringer_area = self.stringer.area
111 Added: cap_area = self.spar.cap_area
112 Added:
113 Added: caps_x = [value for spar in self.spar.x for value in spar]
114 Added: caps_z = [value for spar in self.spar.z for value in spar]
115 Added: stringers_x = self.stringer.x
116 Added: stringers_z = self.stringer.z
117 Added:
118 Added: denominator = float(
119 Added: len(caps_x) * cap_area + len(stringers_x) * stringer_area)
120 Added:
121 Added: centroid_x = float(
122 Added: sum([x * cap_area for x in caps_x]) +
123 Added: sum([x * stringer_area for x in stringers_x]))
124 Added: centroid_x = centroid_x / denominator
125 Added:
126 Added: centroid_z = float(
127 Added: sum([z * cap_area for z in caps_z]) +
128 Added: sum([z * stringer_area for z in stringers_z]))
129 Added: centroid_z = centroid_z / denominator
130 Added:
131 Added: return (centroid_x, centroid_z)
132 Added:
133 Added: def get_inertia_terms(self):
134 Added: """Obtain all inertia terms."""
135 Added: stringer_area = self.stringer.area
136 Added: cap_area = self.spar.cap_area
137 Added:
138 Added: # Adds upper and lower components' coordinates to list
139 Added: x_stringers = self.stringer.x
140 Added: z_stringers = self.stringer.z
141 Added: x_spars = self.spar.x[:][0] + self.spar.x[:][1]
142 Added: z_spars = self.spar.z[:][0] + self.spar.z[:][1]
143 Added: stringer_count = range(len(x_stringers))
144 Added: spar_count = range(len(self.spar.x))
145 Added:
146 Added: # I_x is the sum of the contributions of the spar caps and stringers
147 Added: # TODO: replace list indices with dictionary value
148 Added: I_x = sum([
149 Added: cap_area * (z_spars[i] - self.centroid[1])**2 for i in spar_count
150 Added: ])
151 Added: I_x += sum([
152 Added: stringer_area * (z_stringers[i] - self.centroid[1])**2
153 Added: for i in stringer_count
154 Added: ])
155 Added:
156 Added: I_z = sum([
157 Added: cap_area * (x_spars[i] - self.centroid[0])**2 for i in spar_count
158 Added: ])
159 Added: I_z += sum([
160 Added: stringer_area * (x_stringers[i] - self.centroid[0])**2
161 Added: for i in stringer_count
162 Added: ])
163 Added:
164 Added: I_xz = sum([
165 Added: cap_area * (x_spars[i] - self.centroid[0]) *
166 Added: (z_spars[i] - self.centroid[1]) for i in spar_count
167 Added: ])
168 Added: I_xz += sum([
169 Added: stringer_area * (x_stringers[i] - self.centroid[0]) *
170 Added: (z_stringers[i] - self.centroid[1]) for i in stringer_count
171 Added: ])
172 Added: return (I_x, I_z, I_xz)
173 Added:
174 Added: def get_dx(self, component):
175 Added: return [x - self.centroid[0] for x in component.x_start]
176 Added:
177 Added: def get_dz(self, component):
178 Added: return [x - self.centroid[1] for x in component.x_start]
179 Added:
180 Added: def get_dP(self, xDist, zDist, V_x, V_z, area):
181 Added: I_x = self.I_['x']
182 Added: I_z = self.I_['z']
183 Added: I_xz = self.I_['xz']
184 Added: denom = float(I_x * I_z - I_xz**2)
185 Added: z = float()
186 Added: for _ in range(len(xDist)):
187 Added: z += float(-area * xDist[_] * (I_x * V_x - I_xz * V_z) / denom -
188 Added: area * zDist[_] * (I_z * V_z - I_xz * V_x) / denom)
189 Added: return z
190 Added:
191 Added: def analysis(self, V_x, V_z):
192 Added: """Perform all analysis calculations and store in class instance."""
193 Added: self.drag = self.get_drag(10)
194 Added: self.lift_rectangular = self.get_lift_rectangular(13.7)
195 Added: self.lift_elliptical = self.get_lift_elliptical(15)
196 Added: self.lift_total = self.get_lift_total()
197 Added: self.mass_dist = self.get_mass_distribution(self.mass_total)
198 Added: self.centroid = self.get_centroid()
199 Added: self.I_['x'] = self.get_inertia_terms()[0]
200 Added: self.I_['z'] = self.get_inertia_terms()[1]
201 Added: self.I_['xz'] = self.get_inertia_terms()[2]
202 Added: spar_dx = self.get_dx(self.spar)
203 Added: spar_dz = self.get_dz(self.spar)
204 Added: self.spar.dP_x = self.get_dP(spar_dx, spar_dz, V_x, 0,
205 Added: self.spar.cap_area)
206 Added: self.spar.dP_z = self.get_dP(spar_dx, spar_dz, 0, V_z,
207 Added: self.spar.cap_area)
208 Added: return None
209 Added:
210 Added:
211 Added: def plot_geom(evaluator):
212 Added: """This function plots analysis results over the airfoil's geometry."""
213 Added: # Plot chord
214 Added: x_chord = [0, evaluator.chord]
215 Added: y_chord = [0, 0]
216 Added: plt.plot(x_chord, y_chord, linewidth='1')
217 Added: # Plot quarter chord
218 Added: plt.plot(evaluator.chord / 4,
219 Added: 0,
220 Added: '.',
221 Added: color='g',
222 Added: markersize=24,
223 Added: label='Quarter-chord')
224 Added: # Plot airfoil surfaces
225 Added: x = [0.98 * x for x in evaluator.airfoil.x]
226 Added: y = [0.98 * z for z in evaluator.airfoil.z]
227 Added: plt.fill(x, y, color='w', linewidth='1', fill=False)
228 Added: x = [1.02 * x for x in evaluator.airfoil.x]
229 Added: y = [1.02 * z for z in evaluator.airfoil.z]
230 Added: plt.fill(x, y, color='b', linewidth='1', fill=False)
231 Added:
232 Added: # Plot spars
233 Added: try:
234 Added: for _ in range(len(evaluator.spar.x)):
235 Added: x = (evaluator.spar.x[_])
236 Added: y = (evaluator.spar.z[_])
237 Added: plt.plot(x, y, '-', color='b')
238 Added: except AttributeError:
239 Added: print('No spars to plot.')
240 Added: # Plot stringers
241 Added: try:
242 Added: for _ in range(0, len(evaluator.stringer.x)):
243 Added: x = evaluator.stringer.x[_]
244 Added: y = evaluator.stringer.z[_]
245 Added: plt.plot(x, y, '.', color='y', markersize=12)
246 Added: except AttributeError:
247 Added: print('No stringers to plot.')
248 Added:
249 Added: # Plot centroid
250 Added: x = evaluator.centroid[0]
251 Added: y = evaluator.centroid[1]
252 Added: plt.plot(x, y, '.', color='r', markersize=24, label='centroid')
253 Added:
254 Added: # Graph formatting
255 Added: plt.xlabel('X axis')
256 Added: plt.ylabel('Z axis')
257 Added:
258 Added: plot_bound = max(evaluator.airfoil.x)
259 Added: plt.xlim(-0.10 * plot_bound, 1.10 * plot_bound)
260 Added: plt.ylim(-(1.10 * plot_bound / 2), (1.10 * plot_bound / 2))
261 Added: plt.gca().set_aspect('equal', adjustable='box')
262 Added: plt.gca().legend()
263 Added: plt.grid(axis='both', linestyle=':', linewidth=1)
264 Added: plt.show()
265 Added: return None
266 Added:
267 Added:
268 Added: def plot_lift(evaluator):
269 Added: x = range(evaluator.semi_span)
270 Added: y_1 = evaluator.lift_rectangular
271 Added: y_2 = evaluator.lift_elliptical
272 Added: y_3 = evaluator.lift_total
273 Added: plt.plot(x, y_1, '.', color='b', markersize=4, label='Rectangular lift')
274 Added: plt.plot(x, y_2, '.', color='g', markersize=4, label='Elliptical lift')
275 Added: plt.plot(x, y_3, '.', color='r', markersize=4, label='Total lift')
276 Added:
277 Added: # Graph formatting
278 Added: plt.xlabel('Semi-span location')
279 Added: plt.ylabel('Lift')
280 Added:
281 Added: plt.gca().legend()
282 Added: plt.grid(axis='both', linestyle=':', linewidth=1)
283 Added: plt.show()
284 Added: return None
example_airfoil.py
index 00000000..292c1e99 000000..100644
@@ -0,0 +1,82 @@
1 Added: """This example illustrates the usage of creator, evaluator and generator.
2 Added:
3 Added: All the steps of airfoil creation & evaluation are detailed here;
4 Added: furthermore, the generator.py module contains certain presets
5 Added: (default airfoils).
6 Added:
7 Added: Create an airfoil;
8 Added: Evaluate an airfoil;
9 Added: Generate a population of airfoils & optimize.
10 Added: """
11 Added:
12 Added: from resources import materials as mt
13 Added: from creator import wing, fuselage, propulsion
14 Added: # from evaluator import
15 Added: # from generator import
16 Added:
17 Added: import time
18 Added: start_time = time.time()
19 Added:
20 Added: # Airfoil dimensions (in)
21 Added: NACA_NUM = 2412
22 Added:
23 Added: # Thicknesses
24 Added: SPAR_THICKNESS = 0.4
25 Added: SKIN_THICKNESS = 0.1
26 Added:
27 Added: # Component masses (lbs)
28 Added: AIRFOIL_MASS = 10
29 Added: SPAR_MASS = 10
30 Added: STRINGER_MASS = 5
31 Added:
32 Added: # Area (sqin)
33 Added: SPAR_CAP_AREA = 0.3
34 Added: STRINGER_AREA = 0.1
35 Added:
36 Added: # Amount of stringers
37 Added: TOP_STRINGERS = 6
38 Added: BOTTOM_STRINGERS = 4
39 Added: NOSE_TOP_STRINGERS = 3
40 Added: NOSE_BOTTOM_STRINGERS = 5
41 Added:
42 Added: SAVE_PATH = '/home/blendux/Projects/Aircraft_Studio/save'
43 Added:
44 Added: # Create airfoil instance
45 Added: af = wing.Airfoil(68, 150, mt.aluminium)
46 Added: af.add_naca(NACA_NUM)
47 Added: # af.info_print(2)
48 Added: af.info_save(SAVE_PATH, 'foo_name')
49 Added:
50 Added: # Create spar instances
51 Added: af.spar1 = wing.Spar(af, 0.23, mt.aluminium)
52 Added: af.spar2 = wing.Spar(af, 0.57, mt.aluminium)
53 Added: # af.spar1.info_print(2)
54 Added: # af.spar2.info_print(2)
55 Added: af.spar1.info_save(SAVE_PATH, 'spar1')
56 Added: af.spar2.info_save(SAVE_PATH, 'spar2')
57 Added:
58 Added: # # Create stringer instance
59 Added: # af.stringer = wing.Stringer()
60 Added: # # Compute the stringer coordinates from their quantity in each zone
61 Added: # af.stringer.add_coord(af, [af.spar1, af.spar2], NOSE_TOP_STRINGERS, TOP_STRINGERS,
62 Added: # NOSE_BOTTOM_STRINGERS, BOTTOM_STRINGERS)
63 Added: # af.stringer.add_area(STRINGER_AREA)
64 Added: # af.stringer.add_webs(SKIN_THICKNESS)
65 Added: # af.stringer.info_print(2)
66 Added: # af.stringer.info_save(SAVE_PATH, 'foo_name')
67 Added:
68 Added: # Plot components with matplotlib
69 Added: wing.plot_geom(af, [af.spar1, af.spar2], None)
70 Added:
71 Added: # Evaluator object contains airfoil analysis results.
72 Added: # eval = evaluator.Evaluator(af)
73 Added: # The analysis is performed in the evaluator.py module.
74 Added: # eval.analysis(1, 1)
75 Added: # eval.info_print(2)
76 Added: # eval.info_save(SAVE_PATH, 'foo_name')
77 Added: # evaluator.plot_geom(eval)
78 Added: # evaluator.plot_lift(eval)
79 Added:
80 Added: # Final execution time
81 Added: final_time = time.time() - start_time
82 Added: print(f"--- {round(final_time, 4)}s seconds ---")
generator.py
index 00000000..02138283 000000..100644
@@ -0,0 +1,64 @@
1 Added: # This file is part of Marius Peter's airfoil analysis package (this program).
2 Added: #
3 Added: # This program is free software: you can redistribute it and/or modify
4 Added: # it under the terms of the GNU General Public License as published by
5 Added: # the Free Software Foundation, either version 3 of the License, or
6 Added: # (at your option) any later version.
7 Added: #
8 Added: # This program is distributed in the hope that it will be useful,
9 Added: # but WITHOUT ANY WARRANTY; without even the implied warranty of
10 Added: # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 Added: # GNU General Public License for more details.
12 Added: #
13 Added: # You should have received a copy of the GNU General Public License
14 Added: # along with this program. If not, see <https://www.gnu.org/licenses/>.
15 Added: """
16 Added: The generator.py module contains a single Population class,
17 Added: which represents a collection of randomized airfoils.
18 Added: """
19 Added:
20 Added: from tools import creator
21 Added:
22 Added:
23 Added: def default_airfoil():
24 Added: """Generate the default airfoil."""
25 Added: airfoil = creator.Airfoil.from_dimensions(100, 200)
26 Added: airfoil.add_naca(2412)
27 Added: airfoil.add_mass(10)
28 Added:
29 Added: airfoil.spar = creator.Spar()
30 Added: airfoil.spar.add_coord(airfoil, 0.23)
31 Added: airfoil.spar.add_coord(airfoil, 0.57)
32 Added: airfoil.spar.add_spar_caps(0.3)
33 Added: airfoil.spar.add_mass(10)
34 Added: airfoil.spar.add_webs(0.4)
35 Added:
36 Added: airfoil.stringer = creator.Stringer()
37 Added: airfoil.stringer.add_coord(airfoil, 3, 6, 5, 4)
38 Added: airfoil.stringer.add_area(0.1)
39 Added: airfoil.stringer.add_mass(5)
40 Added: airfoil.stringer.add_webs(0.1)
41 Added:
42 Added: return airfoil
43 Added:
44 Added:
45 Added: class Population(creator.Airfoil):
46 Added: """Collection of random airfoils."""
47 Added:
48 Added: def __init__(self, size):
49 Added: af = creator.Airfoil
50 Added: # print(af)
51 Added: self.size = size
52 Added: self.gen_number = 0 # incremented for every generation
53 Added:
54 Added: def mutate(self, prob_mt):
55 Added: """Randomly mutate the genes of prob_mt % of the population."""
56 Added:
57 Added: def crossover(self, prob_cx):
58 Added: """Combine the genes of prob_cx % of the population."""
59 Added:
60 Added: def reproduce(self, prob_rp):
61 Added: """Pass on the genes of the fittest prob_rp % of the population."""
62 Added:
63 Added: def fitness():
64 Added: """Rate the fitness of an individual on a relative scale (0-100)"""
gui.py
index 00000000..2eea281a 000000..100644
@@ -0,0 +1,81 @@
1 Added: from tools import creator, evaluator, generator
2 Added: # import creator
3 Added: # import evaluator
4 Added: # import generator
5 Added: import tkinter as tk
6 Added: import tkinter.ttk as ttk
7 Added:
8 Added: from matplotlib.backends.backend_tkagg import (
9 Added: FigureCanvasTkAgg, NavigationToolbar2Tk)
10 Added:
11 Added:
12 Added: class MainWindow(tk.Frame):
13 Added: """Main editor window."""
14 Added:
15 Added: def __init__(self, *args, **kwargs):
16 Added: tk.Frame.__init__(self, *args, **kwargs)
17 Added: root = tk.Tk()
18 Added: root.wm_title('MAE 154B - Airfoil Design, Evaluation, Optimization')
19 Added:
20 Added: # self.button = tk.Button(self, text="Create new window",
21 Added: # command=self.create_window)
22 Added: # self.button.pack(side="top")
23 Added: frame_1 = ttk.Frame(root)
24 Added: l_naca, e_naca = new_field(frame_1, 'naca')
25 Added: l_chord, e_chord = new_field(frame_1, 'chord')
26 Added: l_semi_span, e_semi_span = new_field(frame_1, 'semi_span')
27 Added: af = generator.default_airfoil()
28 Added: # Graph window
29 Added: frame_2 = ttk.Frame(root)
30 Added: fig, ax = creator.plot_geom(af, False)
31 Added: plot = FigureCanvasTkAgg(fig, frame_2)
32 Added: # plot.draw()
33 Added: toolbar = NavigationToolbar2Tk(plot, frame_2)
34 Added: # toolbar.update()
35 Added:
36 Added: l_naca.grid(row=0, column=0)
37 Added: e_naca.grid(row=0, column=1, padx=4)
38 Added: # b_naca.grid(row=0, column=2)
39 Added: l_chord.grid(row=1, column=0)
40 Added: e_chord.grid(row=1, column=1, padx=4)
41 Added: l_semi_span.grid(row=2, column=0, padx=4)
42 Added: e_semi_span.grid(row=2, column=1, padx=4)
43 Added: frame_1.pack(side=tk.LEFT)
44 Added: # Graph window
45 Added: plot.get_tk_widget().pack(expand=1, fill=tk.BOTH)
46 Added: toolbar.pack()
47 Added: frame_2.pack(side=tk.LEFT)
48 Added:
49 Added: def create_window(self):
50 Added: self.counter += 1
51 Added: window = tk.Toplevel(self)
52 Added: window.wm_title("Window #%s" % self.counter)
53 Added: label = tk.Label(window, text="This is window #%s" % self.counter)
54 Added: label.pack(side="top", fill="both", expand=True, padx=100, pady=100)
55 Added:
56 Added:
57 Added: def new_field(parent, name):
58 Added: """Add a new user input field."""
59 Added:
60 Added: label = ttk.Label(parent, text=name)
61 Added: entry = ttk.Entry(parent)
62 Added: return label, entry
63 Added:
64 Added:
65 Added: def set_naca(name):
66 Added: naca_num = name.get()
67 Added: print(naca_num)
68 Added:
69 Added:
70 Added: def set_chord(name):
71 Added: chord = name.get()
72 Added: print(chord)
73 Added:
74 Added:
75 Added: def set_semi_span(name):
76 Added: semi_span = name.get()
77 Added: print(semi_span)
78 Added:
79 Added:
80 Added: # plot.get_tk_widget().pack()
81 Added: MainWindow().mainloop()
resources/materials.py
index 00000000..480b518d 000000..100644
@@ -0,0 +1,8 @@
1 Added: aluminium = {
2 Added: "name": "aluminium",
3 Added: "category": "metal",
4 Added: "density": 2.70,
5 Added: "mod_young": 70E9,
6 Added: "mod_shear": 26E9,
7 Added: "mod_bulk": 76E9
8 Added: }
wing_scripts/eye_beam_example.m
index 00000000..70b4d921 000000..100644
@@ -0,0 +1,70 @@
1 Added: % Bending/Shear stress example
2 Added: close all;
3 Added:
4 Added: length = 20; % in
5 Added: force = 10000; %lbs
6 Added:
7 Added: %eye-beam dimensions
8 Added:
9 Added: max_width = 4; % in
10 Added: min_width = 1; % in
11 Added: y_max = 4; % in
12 Added: center_y = 2; % in
13 Added:
14 Added:
15 Added: %max bending moment at the root...
16 Added:
17 Added: M = force*length;
18 Added:
19 Added: I = min_width*(2*center_y)^3/12 + 2*( max_width*(y_max-center_y)^3/12 + ...
20 Added: max_width*(y_max-center_y)*((y_max+center_y)/2)^2);
21 Added:
22 Added: sigma_max = M * y_max / I;
23 Added:
24 Added:
25 Added: % solve for shear stress distribution
26 Added: % V / (I * t) * int(y*da)
27 Added:
28 Added: % Point 1: evaluated at location just before thickness changes from 4 to 1 in
29 Added: tempCoeff = force / (I * max_width);
30 Added: int_y_da = ((y_max+center_y)/2) * max_width*(y_max-center_y);
31 Added: shear_1 = tempCoeff*int_y_da;
32 Added:
33 Added: % Point 2: evaluated at location just after thickness changes from 4 to 1 in
34 Added: tempCoeff = force / (I * min_width);
35 Added: shear_2 = tempCoeff*int_y_da;
36 Added:
37 Added:
38 Added: % Point 3: evaluated at center of beam
39 Added: tempCoeff = force / (I * min_width);
40 Added: int_y_da = (center_y/2) * min_width*center_y;
41 Added: shear_3 = shear_2+tempCoeff*int_y_da;
42 Added:
43 Added: %evaluating continous integral for width of 4..
44 Added: int_y_da_4 = force / (I * max_width)*4*(y_max^2/2 - (center_y:.1:y_max).^2/2);
45 Added:
46 Added: %evaluating continous integral for width of 1..
47 Added: int_y_da_1 = shear_2 + force / (I * min_width)*1*(center_y^2/2 - (0:.1:center_y).^2/2);
48 Added:
49 Added: figure; grid on; hold on;set(gcf,'color',[1 1 1]);
50 Added: plot(int_y_da_4,center_y:.1:y_max,'linewidth',2)
51 Added: plot(int_y_da_1,0:.1:center_y,'linewidth',2)
52 Added: plot(int_y_da_1,0:-.1:-center_y,'linewidth',2)
53 Added: plot(int_y_da_4,-center_y:-.1:-y_max,'linewidth',2)
54 Added: plot([shear_1 shear_2],[center_y center_y],'linewidth',2)
55 Added: plot([shear_1 shear_2],[-center_y -center_y],'linewidth',2)
56 Added:
57 Added: plot(shear_1,center_y,'o')
58 Added: plot(shear_2,center_y,'o')
59 Added: plot(shear_3,0,'o')
60 Added: plot(shear_2,-center_y,'o')
61 Added: plot(shear_1,-center_y,'o')
62 Added:
63 Added: xlabel('shear stress (lb/in^2)','fontsize',16,'fontweight','bold');ylabel('Distance from Center (in)','fontsize',16,'fontweight','bold')
64 Added: set(gca,'FontSize',16,'fontweight','bold');
65 Added:
66 Added:
67 Added: figure; grid on; hold on;set(gcf,'color',[1 1 1]);
68 Added: plot([0 4 4 2.5 2.5 4 4 0 0 1.5 1.5 0 0],[4 4 2 2 -2 -2 -4 -4 -2 -2 2 2 4],'linewidth',2)
69 Added:
70 Added:
wing_scripts/get_dp.m
index 00000000..2a3281d4 000000..100644
@@ -0,0 +1,4 @@
1 Added: function z = get_dp(xDist,zDist,Vx,Vz,Ix,Iz,Ixz,A)
2 Added:
3 Added: denom = (Ix*Iz-Ixz^2);
4 Added: z = -A*xDist*(Ix*Vx-Ixz*Vz)/denom - A*zDist*(Iz*Vz-Ixz*Vx)/denom;
wing_scripts/get_ds.m
index 00000000..2f0eb9d5 000000..100644
@@ -0,0 +1,20 @@
1 Added: function ds = get_ds(xi,xf,u)
2 Added:
3 Added: dist = 0;
4 Added: numSteps = 10;
5 Added: dx = (xf-xi)/numSteps;
6 Added: z0 = get_z(xi,u);
7 Added: x0 = xi;
8 Added: for i=1:10
9 Added: tempX = x0+dx;
10 Added: if tempX > 0
11 Added: tempZ = get_z(tempX,u);
12 Added: else
13 Added: tempZ = 0;
14 Added: end
15 Added: dist = dist + (dx^2+(tempZ-z0)^2)^.5;
16 Added: z0 = tempZ;
17 Added: x0 = tempX;
18 Added: end
19 Added:
20 Added: ds =dist;
wing_scripts/get_int.m
index 00000000..edbfda38 000000..100644
@@ -0,0 +1,35 @@
1 Added: function z = get_int(xi,xf,u)
2 Added:
3 Added: M = 0.02;
4 Added: P = 0.4;
5 Added: T = 0.12;
6 Added: a0 = 0.2969;
7 Added: a1 = -0.126;
8 Added: a2 = -0.3516;
9 Added: a3 = 0.2843;
10 Added: a4 = -0.1015;
11 Added:
12 Added:
13 Added: %evaluate the integral of camber line, depending on xi and xf related to P
14 Added:
15 Added: if xf <P
16 Added: intCamb = M/P^2*(2*P*xf^2/2 - xf^3/3) - M/P^2*(2*P*xi^2/2 - xi^3/3);
17 Added: elseif xi<P
18 Added: intCamb = (M/(1-P)^2)*((1 - 2*P)*xf +2*P*xf^2/2 - xf^3/3) - (M/(1-P)^2)*((1 - 2*P)*P +2*P*P^2/2 - P^3/3);
19 Added: intCamb = intCamb + M/P^2*(2*P*P^2/2 - P^3/3) - M/P^2*(2*P*xi^2/2 - xi^3/3);
20 Added: else
21 Added: intCamb = (M/(1-P)^2)*((1 - 2*P)*xf +2*P*xf^2/2 - xf^3/3) - (M/(1-P)^2)*((1 - 2*P)*xi +2*P*xi^2/2 - xi^3/3);
22 Added: end
23 Added:
24 Added: % do integral on thickness line
25 Added: %z_thickness = (T/0.2)*(a0*x^.5+a1*x+a2*x^2+a3*x^3+a4*x^4);
26 Added:
27 Added: intThickness = (T/0.2)*(a0*xf^1.5/1.5 + a1*xf^2/2 + a2*xf^3/3 + a3*xf^4/4 +a4*xf^5/5);
28 Added: intThickness = intThickness - (T/0.2)*(a0*xi^1.5/1.5 + a1*xi^2/2 + a2*xi^3/3 + a3*xi^4/4 +a4*xi^5/5);
29 Added:
30 Added: % combine both integral results to get total integral
31 Added: if u == 1
32 Added: z = intCamb + intThickness;
33 Added: else
34 Added: z = abs(intCamb - intThickness);
35 Added: end
wing_scripts/get_z.m
index 00000000..5387b52d 000000..100644
@@ -0,0 +1,34 @@
1 Added: function z = get_z(x,u)
2 Added:
3 Added:
4 Added:
5 Added: if (x < 0 )
6 Added: disp('invalid X')
7 Added: end
8 Added:
9 Added: M = 0.02;
10 Added: P = 0.4;
11 Added: T = 0.12;
12 Added: a0 = 0.2969;
13 Added: a1 = -0.126;
14 Added: a2 = -0.3516;
15 Added: a3 = 0.2843;
16 Added: a4 = -0.1015;
17 Added:
18 Added: if x <P
19 Added: z_camber = M/P^2*(2*P*x - x^2);
20 Added: else
21 Added: z_camber = (M/(1-P)^2)*(1 - 2*P +2*P*x - x^2);
22 Added: end
23 Added:
24 Added: %z_camber = M/P^2*(2*P*x - x^2);
25 Added: z_thickness = (T/0.2)*(a0*x^.5+a1*x+a2*x^2+a3*x^3+a4*x^4);
26 Added:
27 Added: if u==1
28 Added: z = z_camber + z_thickness;
29 Added: else
30 Added: z = z_camber - z_thickness;
31 Added: end
32 Added:
33 Added:
34 Added:
wing_scripts/my_progress.m
index 00000000..e87ef238 000000..100644
@@ -0,0 +1,459 @@
1 Added: %wing shear flow
2 Added: clear all;
3 Added: close all;
4 Added:
5 Added: Vx = 1; Vz = 1; My = 1; %test loads will be applied individually
6 Added:
7 Added:
8 Added: %Ixz = -Ixz;
9 Added:
10 Added: %define webs
11 Added:
12 Added: %% web cell 1
13 Added:
14 Added: %upper webs
15 Added: numStringers = numTopStringers;
16 Added: stringerGap = upperStringerGap;
17 Added: webThickness = t_upper;
18 Added: tempStringers = topStringers;
19 Added:
20 Added: for i=1:(numStringers+1)
21 Added: web(i).xStart = sparCaps(1).posX + stringerGap*(i-1);
22 Added: web(i).xEnd = sparCaps(1).posX + stringerGap*(i);
23 Added: web(i).thickness = webThickness;
24 Added: web(i).zStart = get_z(web(i).xStart/chord,1)*chord;
25 Added: web(i).zEnd = get_z(web(i).xEnd/chord,1)*chord;
26 Added: if i==1
27 Added: web(i).dp_area = sparCaps(1).area;
28 Added: web(i).dP_X = 0;
29 Added: web(i).dP_Z = 0;
30 Added: web(i).qPrime_X = 0;
31 Added: web(i).qPrime_Z = 0;
32 Added: else
33 Added: web(i).dp_area = tempStringers(i-1).area;
34 Added: dx = web(i).xStart-centroid.posX; dz = web(i).zStart-centroid.posZ;
35 Added: web(i).dP_X = get_dp(dx,dz,Vx,0,Ix,Iz,Ixz,web(i).dp_area); %just Vx
36 Added: web(i).dP_Z = get_dp(dx,dz,0,Vz,Ix,Iz,Ixz,web(i).dp_area); %just Vz
37 Added: web(i).qPrime_X = web(i-1).qPrime_X - web(i).dP_X;
38 Added: web(i).qPrime_Z = web(i-1).qPrime_Z - web(i).dP_Z;
39 Added: end
40 Added: tempInt = get_int(web(i).xStart/chord,web(i).xEnd/chord,1)*chord^2; %integral of airfoil function
41 Added: triangle1 = abs( (web(i).xStart - sparCaps(1).posX)*web(i).zStart/2);
42 Added: triangle2 = abs((web(i).xEnd - sparCaps(1).posX)*web(i).zEnd/2);
43 Added: web(i).Area = tempInt + triangle1 - triangle2;
44 Added: web(i).ds = get_ds(web(i).xStart/chord,web(i).xEnd/chord,1)*chord;
45 Added: web(i).dS_over_t = web(i).ds / web(i).thickness;
46 Added:
47 Added: web(i).q_dS_over_t_X = web(i).qPrime_X * web(i).dS_over_t;
48 Added: web(i).q_dS_over_t_Z = web(i).qPrime_Z * web(i).dS_over_t;
49 Added: web(i).two_A_qprime_X = 2*web(i).Area*web(i).qPrime_X;
50 Added: web(i).two_A_qprime_Z = 2*web(i).Area*web(i).qPrime_Z;
51 Added: web(i).qp_dx_X = web(i).qPrime_X *(web(i).xEnd-web(i).xStart);
52 Added: web(i).qp_dx_Z = web(i).qPrime_Z *(web(i).xEnd-web(i).xStart);
53 Added: web(i).qp_dz_X = web(i).qPrime_X *(web(i).zEnd-web(i).zStart);
54 Added: web(i).qp_dz_Z = web(i).qPrime_Z *(web(i).zEnd-web(i).zStart);
55 Added: end
56 Added: webTop = web;
57 Added: web = [];
58 Added:
59 Added: %rear spar
60 Added: i=1;
61 Added: web(i).xStart = sparCaps(3).posX;
62 Added: web(i).xEnd = sparCaps(4).posX;
63 Added: web(i).thickness = t_rearSpar;
64 Added: web(i).zStart = sparCaps(3).posZ;
65 Added: web(i).zEnd = sparCaps(4).posZ;
66 Added: web(i).dp_area = sparCaps(3).area;
67 Added: dx = web(i).xStart-centroid.posX; dz = web(i).zStart-centroid.posZ;
68 Added: web(i).dP_X = get_dp(dx,dz,Vx,0,Ix,Iz,Ixz,web(i).dp_area);
69 Added: web(i).dP_Z = get_dp(dx,dz,0,Vz,Ix,Iz,Ixz,web(i).dp_area);
70 Added: web(i).qPrime_X = webTop(numTopStringers+1).qPrime_X - web(i).dP_X;
71 Added: web(i).qPrime_Z = webTop(numTopStringers+1).qPrime_Z - web(i).dP_Z;
72 Added:
73 Added: web(i).Area = (sparCaps(3).posX-sparCaps(1).posX)*sparCaps(3).posZ/2 + ...
74 Added: abs((sparCaps(3).posX-sparCaps(1).posX)*sparCaps(4).posZ/2);
75 Added: web(i).ds = abs(sparCaps(3).posZ - sparCaps(4).posZ);
76 Added: web(i).dS_over_t = web(i).ds / web(i).thickness;
77 Added:
78 Added: web(i).q_dS_over_t_X = web(i).qPrime_X * web(i).dS_over_t;
79 Added: web(i).q_dS_over_t_Z = web(i).qPrime_Z * web(i).dS_over_t;
80 Added: web(i).two_A_qprime_X = 2*web(i).Area*web(i).qPrime_X;
81 Added: web(i).two_A_qprime_Z = 2*web(i).Area*web(i).qPrime_Z;
82 Added: web(i).qp_dx_X = web(i).qPrime_X *(web(i).xEnd-web(i).xStart);
83 Added: web(i).qp_dx_Z = web(i).qPrime_Z *(web(i).xEnd-web(i).xStart);
84 Added: web(i).qp_dz_X = web(i).qPrime_X *(web(i).zEnd-web(i).zStart);
85 Added: web(i).qp_dz_Z = web(i).qPrime_Z *(web(i).zEnd-web(i).zStart);
86 Added:
87 Added: webRearSpar = web;
88 Added: web = [];
89 Added:
90 Added:
91 Added: %lower webs
92 Added: numStringers = numBottomStringers;
93 Added: stringerGap = lowerStringerGap;
94 Added: webThickness = t_lower;
95 Added: tempStringers = bottomStringers;
96 Added:
97 Added: for i=1:(numStringers+1)
98 Added: web(i).xStart = sparCaps(4).posX - stringerGap*(i-1);
99 Added: web(i).xEnd = sparCaps(4).posX - stringerGap*(i);
100 Added: web(i).thickness = webThickness;
101 Added: web(i).zStart = get_z(web(i).xStart/chord,0)*chord;
102 Added: web(i).zEnd = get_z(web(i).xEnd/chord,0)*chord;
103 Added: dx = web(i).xStart-centroid.posX; dz = web(i).zStart-centroid.posZ;
104 Added: if i==1
105 Added: web(i).dp_area = sparCaps(4).area;
106 Added: web(i).dP_X = get_dp(dx,dz,Vx,0,Ix,Iz,Ixz,web(i).dp_area);
107 Added: web(i).dP_Z = get_dp(dx,dz,0,Vz,Ix,Iz,Ixz,web(i).dp_area);
108 Added: web(i).qPrime_X = webRearSpar.qPrime_X - web(i).dP_X;
109 Added: web(i).qPrime_Z = webRearSpar.qPrime_Z - web(i).dP_Z;
110 Added: else
111 Added: web(i).dp_area = tempStringers(i-1).area;
112 Added: web(i).dP_X = get_dp(dx,dz, Vx,0,Ix,Iz,Ixz,web(i).dp_area);
113 Added: web(i).dP_Z = get_dp(dx,dz, 0,Vz,Ix,Iz,Ixz,web(i).dp_area);
114 Added: web(i).qPrime_X = web(i-1).qPrime_X - web(i).dP_X;
115 Added: web(i).qPrime_Z = web(i-1).qPrime_Z - web(i).dP_Z;
116 Added: end
117 Added:
118 Added: tempInt = get_int(web(i).xEnd/chord,web(i).xStart/chord,0)*chord^2; %integral of airfoil function
119 Added: triangle2 = abs((web(i).xStart - sparCaps(1).posX)*web(i).zStart/2);
120 Added: triangle1 = abs((web(i).xEnd - sparCaps(1).posX)*web(i).zEnd/2);
121 Added: web(i).Area = tempInt + triangle1 - triangle2;
122 Added: web(i).ds = get_ds(web(i).xStart/chord,web(i).xEnd/chord,0)*chord;
123 Added: web(i).dS_over_t = web(i).ds / web(i).thickness;
124 Added:
125 Added: web(i).q_dS_over_t_X = web(i).qPrime_X * web(i).dS_over_t;
126 Added: web(i).q_dS_over_t_Z = web(i).qPrime_Z * web(i).dS_over_t;
127 Added: web(i).two_A_qprime_X = 2*web(i).Area*web(i).qPrime_X;
128 Added: web(i).two_A_qprime_Z = 2*web(i).Area*web(i).qPrime_Z;
129 Added: web(i).qp_dx_X = web(i).qPrime_X*(web(i).xEnd-web(i).xStart);
130 Added: web(i).qp_dx_Z = web(i).qPrime_Z*(web(i).xEnd-web(i).xStart);
131 Added: web(i).qp_dz_X = web(i).qPrime_X*(web(i).zEnd-web(i).zStart);
132 Added: web(i).qp_dz_Z = web(i).qPrime_Z*(web(i).zEnd-web(i).zStart);
133 Added:
134 Added: %web(i).radCurv = ... Example: get_curve(web(i).xStart,web(i).xEnd,1)
135 Added: end
136 Added: webBottom = web;
137 Added: web = [];
138 Added:
139 Added: %front Spar
140 Added: i=1;
141 Added: web(i).xStart = sparCaps(2).posX;
142 Added: web(i).xEnd = sparCaps(1).posX;
143 Added: web(i).thickness = t_frontSpar;
144 Added: web(i).zStart = sparCaps(2).posZ;
145 Added: web(i).zEnd = sparCaps(1).posZ;
146 Added: web(i).dp_area = sparCaps(2).area;
147 Added: dx = web(i).xStart-centroid.posX; dz = web(i).zStart-centroid.posZ;
148 Added: web(i).dP_X = get_dp(dx,dz,Vx,0,Ix,Iz,Ixz,web(i).dp_area);
149 Added: web(i).dP_Z = get_dp(dx,dz,0,Vz,Ix,Iz,Ixz,web(i).dp_area);
150 Added: web(i).qPrime_X = webBottom(numBottomStringers+1).qPrime_X - web(i).dP_X;
151 Added: web(i).qPrime_Z = webBottom(numBottomStringers+1).qPrime_Z - web(i).dP_Z;
152 Added: web(i).Area = 0;
153 Added: web(i).ds = abs(sparCaps(2).posZ - sparCaps(1).posZ);
154 Added: web(i).dS_over_t = web(i).ds / web(i).thickness;
155 Added:
156 Added: web(i).q_dS_over_t_X = web(i).qPrime_X * web(i).dS_over_t;
157 Added: web(i).q_dS_over_t_Z = web(i).qPrime_Z * web(i).dS_over_t;
158 Added: web(i).two_A_qprime_X = 2*web(i).Area*web(i).qPrime_X;
159 Added: web(i).two_A_qprime_Z = 2*web(i).Area*web(i).qPrime_Z;
160 Added: web(i).qp_dx_X = web(i).qPrime_X *(web(i).xEnd-web(i).xStart);
161 Added: web(i).qp_dx_Z = web(i).qPrime_Z *(web(i).xEnd-web(i).xStart);
162 Added: web(i).qp_dz_X = web(i).qPrime_X *(web(i).zEnd-web(i).zStart);
163 Added: web(i).qp_dz_Z = web(i).qPrime_Z *(web(i).zEnd-web(i).zStart);
164 Added:
165 Added: webFrontSpar = web;
166 Added: web = [];
167 Added:
168 Added:
169 Added:
170 Added:
171 Added: %% web cell 2
172 Added:
173 Added: %lower nose webs
174 Added: numStringers = numNoseBottomStringers;
175 Added: stringerGap = lowerNoseStringerGap;
176 Added: webThickness = t_lower_front;
177 Added: tempStringers = noseBottomStringers;
178 Added:
179 Added: for i=1:(numStringers+1)
180 Added: web(i).xStart = sparCaps(2).posX - stringerGap*(i-1);
181 Added: web(i).xEnd = sparCaps(2).posX - stringerGap*(i);
182 Added: web(i).thickness = webThickness;
183 Added: web(i).zStart = get_z(web(i).xStart/chord,0)*chord;
184 Added: web(i).zEnd = get_z(web(i).xEnd/chord,0)*chord;
185 Added: dx = web(i).xStart-centroid.posX; dz = web(i).zStart-centroid.posZ;
186 Added:
187 Added: if i==1
188 Added: web(i).dp_area = sparCaps(2).area;
189 Added: web(i).dP_X = 0;
190 Added: web(i).dP_Z = 0;
191 Added: web(i).qPrime_X = 0;
192 Added: web(i).qPrime_Z = 0;
193 Added: else
194 Added: web(i).dp_area = tempStringers(i-1).area;
195 Added: web(i).dP_X = get_dp(dx,dz,Vx,0,Ix,Iz,Ixz,web(i).dp_area);
196 Added: web(i).dP_Z = get_dp(dx,dz,0,Vz,Ix,Iz,Ixz,web(i).dp_area);
197 Added: web(i).qPrime_X = web(i-1).qPrime_X - web(i).dP_X;
198 Added: web(i).qPrime_Z = web(i-1).qPrime_Z - web(i).dP_Z;
199 Added: end
200 Added: tempInt = get_int(web(i).xEnd/chord,web(i).xStart/chord,0)*chord^2; %integral of airfoil function
201 Added: triangle1 = abs((web(i).xStart - sparCaps(2).posX)*web(i).zStart/2);
202 Added: triangle2 = abs((web(i).xEnd - sparCaps(2).posX)*web(i).zEnd/2);
203 Added: web(i).Area = tempInt + triangle1 - triangle2;
204 Added: web(i).ds = get_ds(web(i).xStart/chord,web(i).xEnd/chord,0)*chord;
205 Added: web(i).dS_over_t = web(i).ds / web(i).thickness;
206 Added:
207 Added: web(i).q_dS_over_t_X = web(i).qPrime_X * web(i).dS_over_t;
208 Added: web(i).q_dS_over_t_Z = web(i).qPrime_Z * web(i).dS_over_t;
209 Added: web(i).two_A_qprime_X = 2*web(i).Area*web(i).qPrime_X;
210 Added: web(i).two_A_qprime_Z = 2*web(i).Area*web(i).qPrime_Z;
211 Added: web(i).qp_dx_X = web(i).qPrime_X *(web(i).xEnd-web(i).xStart);
212 Added: web(i).qp_dx_Z = web(i).qPrime_Z *(web(i).xEnd-web(i).xStart);
213 Added: web(i).qp_dz_X = web(i).qPrime_X *(web(i).zEnd-web(i).zStart);
214 Added: web(i).qp_dz_Z = web(i).qPrime_Z *(web(i).zEnd-web(i).zStart);
215 Added:
216 Added: %web(i).radCurv = ... Example: get_curve(web(i).xStart,web(i).xEnd,1)
217 Added: end
218 Added: webLowerNose = web;
219 Added: web = [];
220 Added:
221 Added: %upper nose webs
222 Added: numStringers = numNoseTopStringers;
223 Added: stringerGap = upperNoseStringerGap;
224 Added: webThickness = t_upper_front;
225 Added: tempStringers = noseTopStringers;
226 Added:
227 Added: for i=1:(numStringers+1)
228 Added: web(i).xStart = stringerGap*(i-1);
229 Added: web(i).xEnd = stringerGap*(i);
230 Added: web(i).thickness = webThickness;
231 Added: web(i).zStart = get_z(web(i).xStart/chord,1)*chord;
232 Added: web(i).zEnd = get_z(web(i).xEnd/chord,1)*chord;
233 Added: dx = web(i).xStart-centroid.posX; dz = web(i).zStart-centroid.posZ;
234 Added: if i==1
235 Added: web(i).dp_area = 0;
236 Added: web(i).dP_X = 0;
237 Added: web(i).dP_Z = 0;
238 Added: web(i).qPrime_X = webLowerNose(numNoseBottomStringers+1).qPrime_X - web(i).dP_X;
239 Added: web(i).qPrime_Z = webLowerNose(numNoseBottomStringers+1).qPrime_Z - web(i).dP_Z;
240 Added: else
241 Added: web(i).dp_area = tempStringers(i-1).area;
242 Added: web(i).dP_X = get_dp(dx,dz,Vx,0,Ix,Iz,Ixz,web(i).dp_area);
243 Added: web(i).dP_Z = get_dp(dx,dz,0,Vz,Ix,Iz,Ixz,web(i).dp_area);
244 Added: web(i).qPrime_X = web(i-1).qPrime_X - web(i).dP_X;
245 Added: web(i).qPrime_Z = web(i-1).qPrime_Z - web(i).dP_Z;
246 Added: end
247 Added: tempInt = get_int(web(i).xStart/chord,web(i).xEnd/chord,1)*chord^2; %integral of airfoil function
248 Added: triangle2 = abs((web(i).xStart - sparCaps(2).posX)*web(i).zStart/2);
249 Added: triangle1 = abs((web(i).xEnd - sparCaps(2).posX)*web(i).zEnd/2);
250 Added: web(i).Area = tempInt + triangle1 - triangle2;
251 Added: web(i).ds = get_ds(web(i).xStart/chord,web(i).xEnd/chord,1)*chord;
252 Added: web(i).dS_over_t = web(i).ds / web(i).thickness;
253 Added:
254 Added: web(i).q_dS_over_t_X = web(i).qPrime_X * web(i).dS_over_t;
255 Added: web(i).q_dS_over_t_Z = web(i).qPrime_Z * web(i).dS_over_t;
256 Added: web(i).two_A_qprime_X = 2*web(i).Area*web(i).qPrime_X;
257 Added: web(i).two_A_qprime_Z = 2*web(i).Area*web(i).qPrime_Z;
258 Added: web(i).qp_dx_X = web(i).qPrime_X *(web(i).xEnd-web(i).xStart);
259 Added: web(i).qp_dx_Z = web(i).qPrime_Z *(web(i).xEnd-web(i).xStart);
260 Added: web(i).qp_dz_X = web(i).qPrime_X *(web(i).zEnd-web(i).zStart);
261 Added: web(i).qp_dz_Z = web(i).qPrime_Z *(web(i).zEnd-web(i).zStart);
262 Added:
263 Added: end
264 Added: webUpperNose = web;
265 Added: web = [];
266 Added:
267 Added:
268 Added: %front Spar
269 Added: i=1;
270 Added: web(i).xStart = sparCaps(1).posX;
271 Added: web(i).xEnd = sparCaps(2).posX;
272 Added: web(i).thickness = t_frontSpar;
273 Added: web(i).zStart = sparCaps(1).posZ;
274 Added: web(i).zEnd = sparCaps(2).posZ;
275 Added: web(i).dp_area = sparCaps(1).area;
276 Added: dx = web(i).xStart-centroid.posX; dz = web(i).zStart-centroid.posZ;
277 Added:
278 Added: web(i).dP_X = get_dp(dx,dz,Vx,0,Ix,Iz,Ixz,web(i).dp_area);
279 Added: web(i).dP_Z = get_dp(dx,dz,0,Vz,Ix,Iz,Ixz,web(i).dp_area);
280 Added: web(i).qPrime_X = webUpperNose(numNoseTopStringers+1).qPrime_X - web(i).dP_X;
281 Added: web(i).qPrime_Z = webUpperNose(numNoseTopStringers+1).qPrime_Z - web(i).dP_Z;
282 Added: web(i).Area = 0;
283 Added: web(i).ds = abs(sparCaps(1).posZ - sparCaps(2).posZ);
284 Added: web(i).dS_over_t = web(i).ds / web(i).thickness;
285 Added: web(i).q_dS_over_t_X = web(i).qPrime_X * web(i).dS_over_t;
286 Added: web(i).q_dS_over_t_Z = web(i).qPrime_Z * web(i).dS_over_t;
287 Added: web(i).two_A_qprime_X = 2*web(i).Area*web(i).qPrime_X;
288 Added: web(i).two_A_qprime_Z = 2*web(i).Area*web(i).qPrime_Z;
289 Added: web(i).qp_dx_X = web(i).qPrime_X *(web(i).xEnd-web(i).xStart);
290 Added: web(i).qp_dx_Z = web(i).qPrime_Z *(web(i).xEnd-web(i).xStart);
291 Added: web(i).qp_dz_X = web(i).qPrime_X *(web(i).zEnd-web(i).zStart);
292 Added: web(i).qp_dz_Z = web(i).qPrime_Z *(web(i).zEnd-web(i).zStart);
293 Added:
294 Added: webFrontSparCell2 = web;
295 Added: web = [];
296 Added:
297 Added:
298 Added: %check that q'*dx sums up to Vx
299 Added:
300 Added: Fx = sum([webTop.qp_dx_X])+webRearSpar.qp_dx_X+ sum([webBottom.qp_dx_X])+webFrontSpar.qp_dx_X; %cell 1
301 Added: Fx = Fx + sum([webLowerNose.qp_dx_X])+ sum([webUpperNose.qp_dx_X]); %cell 2
302 Added: Fx
303 Added: Fz = sum([webTop.qp_dz_X])+webRearSpar.qp_dz_X+ sum([webBottom.qp_dz_X])+webFrontSpar.qp_dz_X; %cell 1
304 Added: Fz = Fz + sum([webLowerNose.qp_dz_X])+ sum([webUpperNose.qp_dz_X]); %cell 2
305 Added: Fz
306 Added:
307 Added: %check that q'*dz sums up to Vz
308 Added:
309 Added:
310 Added: Fx = sum([webTop.qp_dx_Z])+webRearSpar.qp_dx_Z+ sum([webBottom.qp_dx_Z])+webFrontSpar.qp_dx_Z; %cell 1
311 Added: Fx = Fx + sum([webLowerNose.qp_dx_Z])+ sum([webUpperNose.qp_dx_Z]); %cell 2
312 Added: Fx
313 Added: Fz = sum([webTop.qp_dz_Z])+webRearSpar.qp_dz_Z+ sum([webBottom.qp_dz_Z])+webFrontSpar.qp_dz_Z; %cell 1
314 Added: Fz = Fz + sum([webLowerNose.qp_dz_Z])+ sum([webUpperNose.qp_dz_Z]); %cell 2
315 Added: Fz
316 Added:
317 Added: %%
318 Added:
319 Added: % sum up the ds/t and q*ds/t to solve 2 equations, 2 unknowns
320 Added:
321 Added: % [A]*[q1s q2s] = B
322 Added:
323 Added: A11 = sum([webTop.dS_over_t])+webRearSpar.dS_over_t+ sum([webBottom.dS_over_t])+webFrontSpar.dS_over_t;
324 Added: A22 = sum([webLowerNose.dS_over_t])+ sum([webUpperNose.dS_over_t])+webFrontSparCell2.dS_over_t;
325 Added: A12 = -webFrontSpar.dS_over_t;
326 Added: A21 = -webFrontSparCell2.dS_over_t;
327 Added:
328 Added: B1_X = sum([webTop.q_dS_over_t_X])+webRearSpar.q_dS_over_t_X+ sum([webBottom.q_dS_over_t_X])+webFrontSpar.q_dS_over_t_X;
329 Added: B2_X = sum([webLowerNose.q_dS_over_t_X])+ sum([webUpperNose.q_dS_over_t_X])+webFrontSparCell2.q_dS_over_t_X;
330 Added: B1_Z = sum([webTop.q_dS_over_t_Z])+webRearSpar.q_dS_over_t_Z+ sum([webBottom.q_dS_over_t_Z])+webFrontSpar.q_dS_over_t_Z;
331 Added: B2_Z = sum([webLowerNose.q_dS_over_t_Z])+ sum([webUpperNose.q_dS_over_t_Z])+webFrontSparCell2.q_dS_over_t_Z;
332 Added:
333 Added: Amat = [A11 A12; A21 A22];
334 Added: Bmat_X = -[B1_X;B2_X];
335 Added: Bmat_Z = -[B1_Z;B2_Z];
336 Added:
337 Added: qs_X = inv(Amat)*Bmat_X;
338 Added: qs_Z = inv(Amat)*Bmat_Z;
339 Added:
340 Added:
341 Added:
342 Added: sum_2_a_q_X = sum([webTop.two_A_qprime_X])+webRearSpar.two_A_qprime_X+ sum([webBottom.two_A_qprime_X]); %cell 1 qprimes
343 Added: sum_2_a_q_X = sum_2_a_q_X + sum([webLowerNose.two_A_qprime_X])+ sum([webUpperNose.two_A_qprime_X]); %cell 2 qprimes
344 Added: sum_2_a_q_X = sum_2_a_q_X + 2*qs_X(1)*(sum([webTop.Area])+webRearSpar.Area+ sum([webBottom.Area]));
345 Added: sum_2_a_q_X = sum_2_a_q_X + 2*qs_X(2)*(sum([webLowerNose.Area])+ sum([webUpperNose.Area]));
346 Added:
347 Added: sum_2_a_q_Z = sum([webTop.two_A_qprime_Z])+webRearSpar.two_A_qprime_Z+ sum([webBottom.two_A_qprime_Z]); %cell 1 qprimes
348 Added: sum_2_a_q_Z = sum_2_a_q_Z + sum([webLowerNose.two_A_qprime_Z])+ sum([webUpperNose.two_A_qprime_Z]); %cell 2 qprimes
349 Added: sum_2_a_q_Z = sum_2_a_q_Z + 2*qs_Z(1)*(sum([webTop.Area])+webRearSpar.Area+ sum([webBottom.Area]));
350 Added: sum_2_a_q_Z = sum_2_a_q_Z + 2*qs_Z(2)*(sum([webLowerNose.Area])+ sum([webUpperNose.Area]));
351 Added:
352 Added: %shear center
353 Added: sc.posX = sum_2_a_q_Z / Vz + frontSpar*chord;
354 Added: sc.posZ = - sum_2_a_q_X / Vx;
355 Added:
356 Added:
357 Added: % now consider the torque representing shifting the load from the quarter
358 Added: % chord to the SC (need to check signs on these moments)
359 Added:
360 Added: torque_Z = Vz*(sc.posX - 0.25*chord);
361 Added: torque_X = -Vx*sc.posZ;
362 Added:
363 Added:
364 Added: Area1 = sum([webTop.Area]) + webRearSpar.Area + sum([webBottom.Area]);
365 Added: %check area
366 Added: Area1_check = get_int(frontSpar,backSpar,1)*chord^2 + get_int(frontSpar,backSpar,0)*chord^2;
367 Added:
368 Added: Area2 = sum([webLowerNose.Area]) + sum([webUpperNose.Area]);
369 Added: Area2_check = get_int(0,frontSpar,1)*chord^2 + get_int(0,frontSpar,0)*chord^2;
370 Added:
371 Added:
372 Added: %for twist equation (see excel spreadsheet example)
373 Added:
374 Added: q1t_over_q2t = (A22/Area2 + webFrontSpar.dS_over_t/Area1)/(A11/Area1 + webFrontSpar.dS_over_t/Area2);
375 Added:
376 Added: q2t = torque_X/(2*Area1*q1t_over_q2t + 2*Area2);
377 Added: q1t = q2t*q1t_over_q2t;
378 Added: qt_X = [q1t;q2t];
379 Added:
380 Added: q2t = torque_Z/(2*Area1*q1t_over_q2t + 2*Area2);
381 Added: q1t = q2t*q1t_over_q2t;
382 Added: qt_Z = [q1t;q2t];
383 Added:
384 Added:
385 Added:
386 Added: % --- - add up all shear flows: qtot = (qPrime + qs) + qt
387 Added:
388 Added:
389 Added:
390 Added:
391 Added: %--- insert force balance to check total shear flows ---
392 Added:
393 Added: % --- --
394 Added:
395 Added:
396 Added: %end
397 Added:
398 Added: sc
399 Added:
400 Added:
401 Added: %plotting airfoil cross-section
402 Added:
403 Added: xChord = 0:.01:1;
404 Added: xChord = xChord*chord;
405 Added: upperSurface = zeros(1,length(xChord));
406 Added: lowerSurface = zeros(1,length(xChord));
407 Added:
408 Added: for i=1:length(xChord)
409 Added: upperSurface(i) = get_z(xChord(i)/chord,1)*chord;
410 Added: lowerSurface(i) = get_z(xChord(i)/chord,0)*chord;
411 Added: end
412 Added:
413 Added: figure; hold on; axis equal; grid on;
414 Added: %plot(xChord,z_camber,'-')
415 Added: plot(xChord,upperSurface,'-k','linewidth',2)
416 Added: plot(xChord,lowerSurface,'-k','linewidth',2)
417 Added: plot([0 1],[0 0],'--k','linewidth',1)
418 Added:
419 Added:
420 Added: for i = 1:length(webTop)
421 Added: vecX = [frontSpar*chord webTop(i).xStart webTop(i).xEnd];
422 Added: vecZ = [0 webTop(i).zStart webTop(i).zEnd];
423 Added: fill(vecX,vecZ,[0.9 0.9 0.9])
424 Added: end
425 Added:
426 Added: for i = 1:length(webBottom)
427 Added: vecX = [frontSpar*chord webBottom(i).xStart webBottom(i).xEnd];
428 Added: vecZ = [0 webBottom(i).zStart webBottom(i).zEnd];
429 Added: fill(vecX,vecZ,[0.9 0.9 0.9])
430 Added: end
431 Added:
432 Added: for i = 1:length(webUpperNose)
433 Added: vecX = [frontSpar*chord webUpperNose(i).xStart webUpperNose(i).xEnd];
434 Added: vecZ = [0 webUpperNose(i).zStart webUpperNose(i).zEnd];
435 Added: fill(vecX,vecZ,[0.7 0.9 1.0])
436 Added: end
437 Added:
438 Added: for i = 1:length(webLowerNose)
439 Added: vecX = [frontSpar*chord webLowerNose(i).xStart webLowerNose(i).xEnd];
440 Added: vecZ = [0 webLowerNose(i).zStart webLowerNose(i).zEnd];
441 Added: fill(vecX,vecZ,[0.7 0.9 1.0])
442 Added: end
443 Added:
444 Added: vecX = [frontSpar*chord sparCaps(3).posX sparCaps(4).posX];
445 Added: vecZ = [0 sparCaps(3).posZ sparCaps(4).posZ];
446 Added: fill(vecX,vecZ,[0.9 0.9 0.9])
447 Added:
448 Added:
449 Added: sparCapSize = 18;
450 Added: stringerSize = 18;
451 Added: plot([sparCaps(1).posX sparCaps(2).posX],[sparCaps(1).posZ sparCaps(2).posZ],'-k','linewidth',2)
452 Added: plot([sparCaps(3).posX sparCaps(4).posX],[sparCaps(3).posZ sparCaps(4).posZ],'-k','linewidth',2)
453 Added: plot([sparCaps.posX],[sparCaps.posZ],'.b','markersize',sparCapSize)
454 Added: plot([topStringers.posX],[topStringers.posZ],'.r','markersize',stringerSize)
455 Added: plot([bottomStringers.posX],[bottomStringers.posZ],'.r','markersize',stringerSize)
456 Added: plot([noseTopStringers.posX],[noseTopStringers.posZ],'.r','markersize',stringerSize)
457 Added: plot([noseBottomStringers.posX],[noseBottomStringers.posZ],'.r','markersize',stringerSize)
458 Added: plot(centroid.posX,centroid.posZ,'.k','markerSize',18)
459 Added: plot(sc.posX,sc.posZ,'.g','markersize',18)
wing_scripts/stringersBeamExample.m
index 00000000..cd3bcb6f 000000..100644
@@ -0,0 +1,47 @@
1 Added: close all;
2 Added: force = 8000; % lbs
3 Added: stringer_A = 0.5; % in^2
4 Added: thickness = 0.04; % in
5 Added:
6 Added: top_stringers_y = 6; % in
7 Added: middle_stringers_y = 2; % in
8 Added:
9 Added: I = 2*stringer_A*top_stringers_y^2 + 2*stringer_A*middle_stringers_y^2;
10 Added:
11 Added: % solve for shear stress distribution. this calc ignores the thickness of
12 Added: % the web between teh stringers (assumes bending taken by stringers)
13 Added: % V / (I * t) * int(y*da)
14 Added:
15 Added: shear_top_web = force / (I*thickness) * top_stringers_y * stringer_A;
16 Added: shear_middle_web = shear_top_web + (force / (I*thickness) * middle_stringers_y * stringer_A);
17 Added:
18 Added: figure; grid on; hold on;set(gcf,'color',[1 1 1]);
19 Added:
20 Added:
21 Added: plot([shear_top_web shear_top_web],[middle_stringers_y top_stringers_y],'linewidth',2);
22 Added: plot([shear_middle_web shear_middle_web],[-middle_stringers_y middle_stringers_y],'linewidth',2);
23 Added: plot([shear_top_web shear_top_web],[-middle_stringers_y -top_stringers_y],'linewidth',2);
24 Added:
25 Added: plot([0 shear_top_web],[top_stringers_y top_stringers_y],'linewidth',2);
26 Added: plot([0 shear_top_web],[-top_stringers_y -top_stringers_y],'linewidth',2);
27 Added: plot([shear_middle_web shear_top_web],[middle_stringers_y middle_stringers_y],'linewidth',2);
28 Added: plot([shear_middle_web shear_top_web],[-middle_stringers_y -middle_stringers_y],'linewidth',2);
29 Added: xlabel('shear stress (lb/in^2)','fontsize',16,'fontweight','bold');ylabel('Distance from Center (in)','fontsize',16,'fontweight','bold')
30 Added: set(gca,'FontSize',16,'fontweight','bold');
31 Added:
32 Added: %Alternate approach.. compute change in bending stress at each stringer to
33 Added: %find the change in shear load
34 Added:
35 Added: %at top stringer
36 Added: d_sigma = force * top_stringers_y / I; %(lbs/in^2)
37 Added: d_force_top = d_sigma * stringer_A;
38 Added:
39 Added: %at middle stringer..
40 Added: d_sigma = force * middle_stringers_y / I; %(lbs/in^2)
41 Added: d_force_middle = d_force_top + d_sigma*stringer_A;
42 Added:
43 Added: %check if load balances
44 Added: check_load = 2*d_force_top*4 + d_force_middle*4;
45 Added:
46 Added:
47 Added:
wing_scripts/wingAnalysis_190422.m
index 00000000..a7d65e20 000000..100644
@@ -0,0 +1,579 @@
1 Added: %wing shear flow
2 Added: clear all;
3 Added: close all;
4 Added:
5 Added:
6 Added:
7 Added:
8 Added: Vx = 1; Vz = 1; My = 1; %test loads will be applied individually
9 Added:
10 Added: %define a few
11 Added: numTopStringers = 6;
12 Added: numBottomStringers = 8;
13 Added: numNoseTopStringers = 4;
14 Added: numNoseBottomStringers = 4;
15 Added:
16 Added: t_upper = 0.02/12;
17 Added: t_lower = 0.02/12;
18 Added: t_upper_front = 0.02/12;
19 Added: t_lower_front = 0.02/12;
20 Added: t_frontSpar = 0.04/12;
21 Added: t_rearSpar = 0.04/12;
22 Added:
23 Added: frontSpar = 0.2;
24 Added: backSpar = 0.7;
25 Added: chord = 5;
26 Added:
27 Added: sparCaps(1).posX = frontSpar*chord;
28 Added: sparCaps(2).posX = frontSpar*chord;
29 Added: sparCaps(3).posX = backSpar*chord;
30 Added: sparCaps(4).posX = backSpar*chord;
31 Added:
32 Added: sparCaps(1).posZ = get_z(frontSpar,1)*chord;
33 Added: sparCaps(2).posZ = get_z(frontSpar,0)*chord;
34 Added: sparCaps(3).posZ = get_z(backSpar,1)*chord;
35 Added: sparCaps(4).posZ = get_z(backSpar,0)*chord;
36 Added:
37 Added: sparCaps(1).area = .1;
38 Added: sparCaps(2).area = .1;
39 Added: sparCaps(3).area = .1;
40 Added: sparCaps(4).area = .1;
41 Added:
42 Added: upperStringerGap = (sparCaps(3).posX - sparCaps(1).posX)/(numTopStringers + 1);
43 Added: lowerStringerGap = (sparCaps(3).posX - sparCaps(1).posX)/(numBottomStringers + 1);
44 Added: upperNoseStringerGap = (sparCaps(1).posX - 0)/(numNoseTopStringers + 1);
45 Added: lowerNoseStringerGap = (sparCaps(1).posX - 0)/(numNoseBottomStringers + 1);
46 Added:
47 Added:
48 Added: %set stringers spaced evenly along X axis betwen Spars
49 Added: %top Stringers
50 Added: for i=1:numTopStringers
51 Added: topStringers(i).posX = sparCaps(1).posX + upperStringerGap*i;
52 Added: topStringers(i).posZ = get_z(topStringers(i).posX/chord,1)*chord;
53 Added: topStringers(i).area = .1;
54 Added: end
55 Added:
56 Added: %bottom Stringers
57 Added: for i=1:numBottomStringers
58 Added: bottomStringers(i).posX = sparCaps(4).posX - lowerStringerGap*i;
59 Added: bottomStringers(i).posZ = get_z(bottomStringers(i).posX/chord,0)*chord;
60 Added: bottomStringers(i).area = .1;
61 Added:
62 Added: end
63 Added:
64 Added: %nose bottom Stringers
65 Added: for i=1:numNoseBottomStringers
66 Added: noseBottomStringers(i).posX = sparCaps(2).posX - lowerNoseStringerGap*i;
67 Added: noseBottomStringers(i).posZ = get_z(noseBottomStringers(i).posX/chord,0)*chord;
68 Added: noseBottomStringers(i).area = .1;
69 Added: end
70 Added:
71 Added: %nose top Stringers
72 Added: for i=1:numNoseTopStringers
73 Added: noseTopStringers(i).posX = upperNoseStringerGap*i;
74 Added: noseTopStringers(i).posZ = get_z(noseTopStringers(i).posX/chord,1)*chord;
75 Added: noseTopStringers(i).area = .1;
76 Added: end
77 Added:
78 Added:
79 Added: centroid.posX = sum([sparCaps.posX].*[sparCaps.area]) + ...
80 Added: sum([topStringers.posX].*[topStringers.area]) + ...
81 Added: sum([bottomStringers.posX].*[bottomStringers.area]) + ...
82 Added: sum([noseTopStringers.posX].*[noseTopStringers.area]) + ...
83 Added: sum([noseBottomStringers.posX].*[noseBottomStringers.area]);
84 Added:
85 Added: centroid.posX = centroid.posX / ( sum([sparCaps.area]) + sum([topStringers.area]) + ...
86 Added: sum([bottomStringers.area]) + sum([noseTopStringers.area]) + sum([noseBottomStringers.area]));
87 Added:
88 Added: centroid.posZ = sum([sparCaps.posZ].*[sparCaps.area]) + ...
89 Added: sum([topStringers.posZ].*[topStringers.area]) + ...
90 Added: sum([bottomStringers.posZ].*[bottomStringers.area]) + ...
91 Added: sum([noseTopStringers.posZ].*[noseTopStringers.area]) + ...
92 Added: sum([noseBottomStringers.posZ].*[noseBottomStringers.area]);
93 Added:
94 Added: centroid.posZ = centroid.posZ / ( sum([sparCaps.area]) + sum([topStringers.area]) + ...
95 Added: sum([bottomStringers.area]) + sum([noseTopStringers.area]) + sum([noseBottomStringers.area]));
96 Added:
97 Added: %summing contributions for inertia terms
98 Added: Ix = 0; Iz = 0; Ixz = 0;
99 Added:
100 Added: for i=1:4 %spar caps
101 Added: Ix = Ix + sparCaps(i).area*(sparCaps(i).posZ-centroid.posZ)^2;
102 Added: Iz = Iz + sparCaps(i).area*(sparCaps(i).posX-centroid.posX)^2;
103 Added: Ixz = Ixz + sparCaps(i).area*(sparCaps(i).posX-centroid.posX)*(sparCaps(i).posZ-centroid.posZ);
104 Added: end
105 Added:
106 Added:
107 Added: for i=1:numTopStringers %top stringers
108 Added: Ix = Ix + topStringers(i).area*(topStringers(i).posZ-centroid.posZ)^2;
109 Added: Iz = Iz + topStringers(i).area*(topStringers(i).posX-centroid.posX)^2;
110 Added: Ixz = Ixz + topStringers(i).area*(topStringers(i).posX-centroid.posX)*(topStringers(i).posZ-centroid.posZ);
111 Added: end
112 Added: for i=1:numBottomStringers %bottom stringers
113 Added: Ix = Ix + bottomStringers(i).area*(bottomStringers(i).posZ-centroid.posZ)^2;
114 Added: Iz = Iz + bottomStringers(i).area*(bottomStringers(i).posX-centroid.posX)^2;
115 Added: Ixz = Ixz + bottomStringers(i).area*(bottomStringers(i).posX-centroid.posX)*(bottomStringers(i).posZ-centroid.posZ);
116 Added: end
117 Added: for i=1:numNoseTopStringers %nose top stringers
118 Added: Ix = Ix + noseTopStringers(i).area*(noseTopStringers(i).posZ-centroid.posZ)^2;
119 Added: Iz = Iz + noseTopStringers(i).area*(noseTopStringers(i).posX-centroid.posX)^2;
120 Added: Ixz = Ixz + noseTopStringers(i).area*(noseTopStringers(i).posX-centroid.posX)*(noseTopStringers(i).posZ-centroid.posZ);
121 Added: end
122 Added: for i=1:numNoseBottomStringers %nose bottom stringers
123 Added: Ix = Ix + noseBottomStringers(i).area*(noseBottomStringers(i).posZ-centroid.posZ)^2;
124 Added: Iz = Iz + noseBottomStringers(i).area*(noseBottomStringers(i).posX-centroid.posX)^2;
125 Added: Ixz = Ixz + noseBottomStringers(i).area*(noseBottomStringers(i).posX-centroid.posX)*(noseBottomStringers(i).posZ-centroid.posZ);
126 Added: end
127 Added:
128 Added: %Ixz = -Ixz;
129 Added:
130 Added: %define webs
131 Added:
132 Added: %% web cell 1
133 Added:
134 Added: %upper webs
135 Added: numStringers = numTopStringers;
136 Added: stringerGap = upperStringerGap;
137 Added: webThickness = t_upper;
138 Added: tempStringers = topStringers;
139 Added:
140 Added: for i=1:(numStringers+1)
141 Added: web(i).xStart = sparCaps(1).posX + stringerGap*(i-1);
142 Added: web(i).xEnd = sparCaps(1).posX + stringerGap*(i);
143 Added: web(i).thickness = webThickness;
144 Added: web(i).zStart = get_z(web(i).xStart/chord,1)*chord;
145 Added: web(i).zEnd = get_z(web(i).xEnd/chord,1)*chord;
146 Added: if i==1
147 Added: web(i).dp_area = sparCaps(1).area;
148 Added: web(i).dP_X = 0;
149 Added: web(i).dP_Z = 0;
150 Added: web(i).qPrime_X = 0;
151 Added: web(i).qPrime_Z = 0;
152 Added: else
153 Added: web(i).dp_area = tempStringers(i-1).area;
154 Added: dx = web(i).xStart-centroid.posX; dz = web(i).zStart-centroid.posZ;
155 Added: web(i).dP_X = get_dp(dx,dz,Vx,0,Ix,Iz,Ixz,web(i).dp_area); %just Vx
156 Added: web(i).dP_Z = get_dp(dx,dz,0,Vz,Ix,Iz,Ixz,web(i).dp_area); %just Vz
157 Added: web(i).qPrime_X = web(i-1).qPrime_X - web(i).dP_X;
158 Added: web(i).qPrime_Z = web(i-1).qPrime_Z - web(i).dP_Z;
159 Added: end
160 Added: tempInt = get_int(web(i).xStart/chord,web(i).xEnd/chord,1)*chord^2; %integral of airfoil function
161 Added: triangle1 = abs( (web(i).xStart - sparCaps(1).posX)*web(i).zStart/2);
162 Added: triangle2 = abs((web(i).xEnd - sparCaps(1).posX)*web(i).zEnd/2);
163 Added: web(i).Area = tempInt + triangle1 - triangle2;
164 Added: web(i).ds = get_ds(web(i).xStart/chord,web(i).xEnd/chord,1)*chord;
165 Added: web(i).dS_over_t = web(i).ds / web(i).thickness;
166 Added:
167 Added: web(i).q_dS_over_t_X = web(i).qPrime_X * web(i).dS_over_t;
168 Added: web(i).q_dS_over_t_Z = web(i).qPrime_Z * web(i).dS_over_t;
169 Added: web(i).two_A_qprime_X = 2*web(i).Area*web(i).qPrime_X;
170 Added: web(i).two_A_qprime_Z = 2*web(i).Area*web(i).qPrime_Z;
171 Added: web(i).qp_dx_X = web(i).qPrime_X *(web(i).xEnd-web(i).xStart);
172 Added: web(i).qp_dx_Z = web(i).qPrime_Z *(web(i).xEnd-web(i).xStart);
173 Added: web(i).qp_dz_X = web(i).qPrime_X *(web(i).zEnd-web(i).zStart);
174 Added: web(i).qp_dz_Z = web(i).qPrime_Z *(web(i).zEnd-web(i).zStart);
175 Added: end
176 Added: webTop = web;
177 Added: web = [];
178 Added:
179 Added: %rear spar
180 Added: i=1;
181 Added: web(i).xStart = sparCaps(3).posX;
182 Added: web(i).xEnd = sparCaps(4).posX;
183 Added: web(i).thickness = t_rearSpar;
184 Added: web(i).zStart = sparCaps(3).posZ;
185 Added: web(i).zEnd = sparCaps(4).posZ;
186 Added: web(i).dp_area = sparCaps(3).area;
187 Added: dx = web(i).xStart-centroid.posX; dz = web(i).zStart-centroid.posZ;
188 Added: web(i).dP_X = get_dp(dx,dz,Vx,0,Ix,Iz,Ixz,web(i).dp_area);
189 Added: web(i).dP_Z = get_dp(dx,dz,0,Vz,Ix,Iz,Ixz,web(i).dp_area);
190 Added: web(i).qPrime_X = webTop(numTopStringers+1).qPrime_X - web(i).dP_X;
191 Added: web(i).qPrime_Z = webTop(numTopStringers+1).qPrime_Z - web(i).dP_Z;
192 Added:
193 Added: web(i).Area = (sparCaps(3).posX-sparCaps(1).posX)*sparCaps(3).posZ/2 + ...
194 Added: abs((sparCaps(3).posX-sparCaps(1).posX)*sparCaps(4).posZ/2);
195 Added: web(i).ds = abs(sparCaps(3).posZ - sparCaps(4).posZ);
196 Added: web(i).dS_over_t = web(i).ds / web(i).thickness;
197 Added:
198 Added: web(i).q_dS_over_t_X = web(i).qPrime_X * web(i).dS_over_t;
199 Added: web(i).q_dS_over_t_Z = web(i).qPrime_Z * web(i).dS_over_t;
200 Added: web(i).two_A_qprime_X = 2*web(i).Area*web(i).qPrime_X;
201 Added: web(i).two_A_qprime_Z = 2*web(i).Area*web(i).qPrime_Z;
202 Added: web(i).qp_dx_X = web(i).qPrime_X *(web(i).xEnd-web(i).xStart);
203 Added: web(i).qp_dx_Z = web(i).qPrime_Z *(web(i).xEnd-web(i).xStart);
204 Added: web(i).qp_dz_X = web(i).qPrime_X *(web(i).zEnd-web(i).zStart);
205 Added: web(i).qp_dz_Z = web(i).qPrime_Z *(web(i).zEnd-web(i).zStart);
206 Added:
207 Added: webRearSpar = web;
208 Added: web = [];
209 Added:
210 Added:
211 Added: %lower webs
212 Added: numStringers = numBottomStringers;
213 Added: stringerGap = lowerStringerGap;
214 Added: webThickness = t_lower;
215 Added: tempStringers = bottomStringers;
216 Added:
217 Added: for i=1:(numStringers+1)
218 Added: web(i).xStart = sparCaps(4).posX - stringerGap*(i-1);
219 Added: web(i).xEnd = sparCaps(4).posX - stringerGap*(i);
220 Added: web(i).thickness = webThickness;
221 Added: web(i).zStart = get_z(web(i).xStart/chord,0)*chord;
222 Added: web(i).zEnd = get_z(web(i).xEnd/chord,0)*chord;
223 Added: dx = web(i).xStart-centroid.posX; dz = web(i).zStart-centroid.posZ;
224 Added: if i==1
225 Added: web(i).dp_area = sparCaps(4).area;
226 Added: web(i).dP_X = get_dp(dx,dz,Vx,0,Ix,Iz,Ixz,web(i).dp_area);
227 Added: web(i).dP_Z = get_dp(dx,dz,0,Vz,Ix,Iz,Ixz,web(i).dp_area);
228 Added: web(i).qPrime_X = webRearSpar.qPrime_X - web(i).dP_X;
229 Added: web(i).qPrime_Z = webRearSpar.qPrime_Z - web(i).dP_Z;
230 Added: else
231 Added: web(i).dp_area = tempStringers(i-1).area;
232 Added: web(i).dP_X = get_dp(dx,dz, Vx,0,Ix,Iz,Ixz,web(i).dp_area);
233 Added: web(i).dP_Z = get_dp(dx,dz, 0,Vz,Ix,Iz,Ixz,web(i).dp_area);
234 Added: web(i).qPrime_X = web(i-1).qPrime_X - web(i).dP_X;
235 Added: web(i).qPrime_Z = web(i-1).qPrime_Z - web(i).dP_Z;
236 Added: end
237 Added:
238 Added: tempInt = get_int(web(i).xEnd/chord,web(i).xStart/chord,0)*chord^2; %integral of airfoil function
239 Added: triangle2 = abs((web(i).xStart - sparCaps(1).posX)*web(i).zStart/2);
240 Added: triangle1 = abs((web(i).xEnd - sparCaps(1).posX)*web(i).zEnd/2);
241 Added: web(i).Area = tempInt + triangle1 - triangle2;
242 Added: web(i).ds = get_ds(web(i).xStart/chord,web(i).xEnd/chord,0)*chord;
243 Added: web(i).dS_over_t = web(i).ds / web(i).thickness;
244 Added:
245 Added: web(i).q_dS_over_t_X = web(i).qPrime_X * web(i).dS_over_t;
246 Added: web(i).q_dS_over_t_Z = web(i).qPrime_Z * web(i).dS_over_t;
247 Added: web(i).two_A_qprime_X = 2*web(i).Area*web(i).qPrime_X;
248 Added: web(i).two_A_qprime_Z = 2*web(i).Area*web(i).qPrime_Z;
249 Added: web(i).qp_dx_X = web(i).qPrime_X*(web(i).xEnd-web(i).xStart);
250 Added: web(i).qp_dx_Z = web(i).qPrime_Z*(web(i).xEnd-web(i).xStart);
251 Added: web(i).qp_dz_X = web(i).qPrime_X*(web(i).zEnd-web(i).zStart);
252 Added: web(i).qp_dz_Z = web(i).qPrime_Z*(web(i).zEnd-web(i).zStart);
253 Added:
254 Added: %web(i).radCurv = ... Example: get_curve(web(i).xStart,web(i).xEnd,1)
255 Added: end
256 Added: webBottom = web;
257 Added: web = [];
258 Added:
259 Added: %front Spar
260 Added: i=1;
261 Added: web(i).xStart = sparCaps(2).posX;
262 Added: web(i).xEnd = sparCaps(1).posX;
263 Added: web(i).thickness = t_frontSpar;
264 Added: web(i).zStart = sparCaps(2).posZ;
265 Added: web(i).zEnd = sparCaps(1).posZ;
266 Added: web(i).dp_area = sparCaps(2).area;
267 Added: dx = web(i).xStart-centroid.posX; dz = web(i).zStart-centroid.posZ;
268 Added: web(i).dP_X = get_dp(dx,dz,Vx,0,Ix,Iz,Ixz,web(i).dp_area);
269 Added: web(i).dP_Z = get_dp(dx,dz,0,Vz,Ix,Iz,Ixz,web(i).dp_area);
270 Added: web(i).qPrime_X = webBottom(numBottomStringers+1).qPrime_X - web(i).dP_X;
271 Added: web(i).qPrime_Z = webBottom(numBottomStringers+1).qPrime_Z - web(i).dP_Z;
272 Added: web(i).Area = 0;
273 Added: web(i).ds = abs(sparCaps(2).posZ - sparCaps(1).posZ);
274 Added: web(i).dS_over_t = web(i).ds / web(i).thickness;
275 Added:
276 Added: web(i).q_dS_over_t_X = web(i).qPrime_X * web(i).dS_over_t;
277 Added: web(i).q_dS_over_t_Z = web(i).qPrime_Z * web(i).dS_over_t;
278 Added: web(i).two_A_qprime_X = 2*web(i).Area*web(i).qPrime_X;
279 Added: web(i).two_A_qprime_Z = 2*web(i).Area*web(i).qPrime_Z;
280 Added: web(i).qp_dx_X = web(i).qPrime_X *(web(i).xEnd-web(i).xStart);
281 Added: web(i).qp_dx_Z = web(i).qPrime_Z *(web(i).xEnd-web(i).xStart);
282 Added: web(i).qp_dz_X = web(i).qPrime_X *(web(i).zEnd-web(i).zStart);
283 Added: web(i).qp_dz_Z = web(i).qPrime_Z *(web(i).zEnd-web(i).zStart);
284 Added:
285 Added: webFrontSpar = web;
286 Added: web = [];
287 Added:
288 Added:
289 Added:
290 Added:
291 Added: %% web cell 2
292 Added:
293 Added: %lower nose webs
294 Added: numStringers = numNoseBottomStringers;
295 Added: stringerGap = lowerNoseStringerGap;
296 Added: webThickness = t_lower_front;
297 Added: tempStringers = noseBottomStringers;
298 Added:
299 Added: for i=1:(numStringers+1)
300 Added: web(i).xStart = sparCaps(2).posX - stringerGap*(i-1);
301 Added: web(i).xEnd = sparCaps(2).posX - stringerGap*(i);
302 Added: web(i).thickness = webThickness;
303 Added: web(i).zStart = get_z(web(i).xStart/chord,0)*chord;
304 Added: web(i).zEnd = get_z(web(i).xEnd/chord,0)*chord;
305 Added: dx = web(i).xStart-centroid.posX; dz = web(i).zStart-centroid.posZ;
306 Added:
307 Added: if i==1
308 Added: web(i).dp_area = sparCaps(2).area;
309 Added: web(i).dP_X = 0;
310 Added: web(i).dP_Z = 0;
311 Added: web(i).qPrime_X = 0;
312 Added: web(i).qPrime_Z = 0;
313 Added: else
314 Added: web(i).dp_area = tempStringers(i-1).area;
315 Added: web(i).dP_X = get_dp(dx,dz,Vx,0,Ix,Iz,Ixz,web(i).dp_area);
316 Added: web(i).dP_Z = get_dp(dx,dz,0,Vz,Ix,Iz,Ixz,web(i).dp_area);
317 Added: web(i).qPrime_X = web(i-1).qPrime_X - web(i).dP_X;
318 Added: web(i).qPrime_Z = web(i-1).qPrime_Z - web(i).dP_Z;
319 Added: end
320 Added: tempInt = get_int(web(i).xEnd/chord,web(i).xStart/chord,0)*chord^2; %integral of airfoil function
321 Added: triangle1 = abs((web(i).xStart - sparCaps(2).posX)*web(i).zStart/2);
322 Added: triangle2 = abs((web(i).xEnd - sparCaps(2).posX)*web(i).zEnd/2);
323 Added: web(i).Area = tempInt + triangle1 - triangle2;
324 Added: web(i).ds = get_ds(web(i).xStart/chord,web(i).xEnd/chord,0)*chord;
325 Added: web(i).dS_over_t = web(i).ds / web(i).thickness;
326 Added:
327 Added: web(i).q_dS_over_t_X = web(i).qPrime_X * web(i).dS_over_t;
328 Added: web(i).q_dS_over_t_Z = web(i).qPrime_Z * web(i).dS_over_t;
329 Added: web(i).two_A_qprime_X = 2*web(i).Area*web(i).qPrime_X;
330 Added: web(i).two_A_qprime_Z = 2*web(i).Area*web(i).qPrime_Z;
331 Added: web(i).qp_dx_X = web(i).qPrime_X *(web(i).xEnd-web(i).xStart);
332 Added: web(i).qp_dx_Z = web(i).qPrime_Z *(web(i).xEnd-web(i).xStart);
333 Added: web(i).qp_dz_X = web(i).qPrime_X *(web(i).zEnd-web(i).zStart);
334 Added: web(i).qp_dz_Z = web(i).qPrime_Z *(web(i).zEnd-web(i).zStart);
335 Added:
336 Added: %web(i).radCurv = ... Example: get_curve(web(i).xStart,web(i).xEnd,1)
337 Added: end
338 Added: webLowerNose = web;
339 Added: web = [];
340 Added:
341 Added: %upper nose webs
342 Added: numStringers = numNoseTopStringers;
343 Added: stringerGap = upperNoseStringerGap;
344 Added: webThickness = t_upper_front;
345 Added: tempStringers = noseTopStringers;
346 Added:
347 Added: for i=1:(numStringers+1)
348 Added: web(i).xStart = stringerGap*(i-1);
349 Added: web(i).xEnd = stringerGap*(i);
350 Added: web(i).thickness = webThickness;
351 Added: web(i).zStart = get_z(web(i).xStart/chord,1)*chord;
352 Added: web(i).zEnd = get_z(web(i).xEnd/chord,1)*chord;
353 Added: dx = web(i).xStart-centroid.posX; dz = web(i).zStart-centroid.posZ;
354 Added: if i==1
355 Added: web(i).dp_area = 0;
356 Added: web(i).dP_X = 0;
357 Added: web(i).dP_Z = 0;
358 Added: web(i).qPrime_X = webLowerNose(numNoseBottomStringers+1).qPrime_X - web(i).dP_X;
359 Added: web(i).qPrime_Z = webLowerNose(numNoseBottomStringers+1).qPrime_Z - web(i).dP_Z;
360 Added: else
361 Added: web(i).dp_area = tempStringers(i-1).area;
362 Added: web(i).dP_X = get_dp(dx,dz,Vx,0,Ix,Iz,Ixz,web(i).dp_area);
363 Added: web(i).dP_Z = get_dp(dx,dz,0,Vz,Ix,Iz,Ixz,web(i).dp_area);
364 Added: web(i).qPrime_X = web(i-1).qPrime_X - web(i).dP_X;
365 Added: web(i).qPrime_Z = web(i-1).qPrime_Z - web(i).dP_Z;
366 Added: end
367 Added: tempInt = get_int(web(i).xStart/chord,web(i).xEnd/chord,1)*chord^2; %integral of airfoil function
368 Added: triangle2 = abs((web(i).xStart - sparCaps(2).posX)*web(i).zStart/2);
369 Added: triangle1 = abs((web(i).xEnd - sparCaps(2).posX)*web(i).zEnd/2);
370 Added: web(i).Area = tempInt + triangle1 - triangle2;
371 Added: web(i).ds = get_ds(web(i).xStart/chord,web(i).xEnd/chord,1)*chord;
372 Added: web(i).dS_over_t = web(i).ds / web(i).thickness;
373 Added:
374 Added: web(i).q_dS_over_t_X = web(i).qPrime_X * web(i).dS_over_t;
375 Added: web(i).q_dS_over_t_Z = web(i).qPrime_Z * web(i).dS_over_t;
376 Added: web(i).two_A_qprime_X = 2*web(i).Area*web(i).qPrime_X;
377 Added: web(i).two_A_qprime_Z = 2*web(i).Area*web(i).qPrime_Z;
378 Added: web(i).qp_dx_X = web(i).qPrime_X *(web(i).xEnd-web(i).xStart);
379 Added: web(i).qp_dx_Z = web(i).qPrime_Z *(web(i).xEnd-web(i).xStart);
380 Added: web(i).qp_dz_X = web(i).qPrime_X *(web(i).zEnd-web(i).zStart);
381 Added: web(i).qp_dz_Z = web(i).qPrime_Z *(web(i).zEnd-web(i).zStart);
382 Added:
383 Added: end
384 Added: webUpperNose = web;
385 Added: web = [];
386 Added:
387 Added:
388 Added: %front Spar
389 Added: i=1;
390 Added: web(i).xStart = sparCaps(1).posX;
391 Added: web(i).xEnd = sparCaps(2).posX;
392 Added: web(i).thickness = t_frontSpar;
393 Added: web(i).zStart = sparCaps(1).posZ;
394 Added: web(i).zEnd = sparCaps(2).posZ;
395 Added: web(i).dp_area = sparCaps(1).area;
396 Added: dx = web(i).xStart-centroid.posX; dz = web(i).zStart-centroid.posZ;
397 Added:
398 Added: web(i).dP_X = get_dp(dx,dz,Vx,0,Ix,Iz,Ixz,web(i).dp_area);
399 Added: web(i).dP_Z = get_dp(dx,dz,0,Vz,Ix,Iz,Ixz,web(i).dp_area);
400 Added: web(i).qPrime_X = webUpperNose(numNoseTopStringers+1).qPrime_X - web(i).dP_X;
401 Added: web(i).qPrime_Z = webUpperNose(numNoseTopStringers+1).qPrime_Z - web(i).dP_Z;
402 Added: web(i).Area = 0;
403 Added: web(i).ds = abs(sparCaps(1).posZ - sparCaps(2).posZ);
404 Added: web(i).dS_over_t = web(i).ds / web(i).thickness;
405 Added: web(i).q_dS_over_t_X = web(i).qPrime_X * web(i).dS_over_t;
406 Added: web(i).q_dS_over_t_Z = web(i).qPrime_Z * web(i).dS_over_t;
407 Added: web(i).two_A_qprime_X = 2*web(i).Area*web(i).qPrime_X;
408 Added: web(i).two_A_qprime_Z = 2*web(i).Area*web(i).qPrime_Z;
409 Added: web(i).qp_dx_X = web(i).qPrime_X *(web(i).xEnd-web(i).xStart);
410 Added: web(i).qp_dx_Z = web(i).qPrime_Z *(web(i).xEnd-web(i).xStart);
411 Added: web(i).qp_dz_X = web(i).qPrime_X *(web(i).zEnd-web(i).zStart);
412 Added: web(i).qp_dz_Z = web(i).qPrime_Z *(web(i).zEnd-web(i).zStart);
413 Added:
414 Added: webFrontSparCell2 = web;
415 Added: web = [];
416 Added:
417 Added:
418 Added: %check that q'*dx sums up to Vx
419 Added:
420 Added: Fx = sum([webTop.qp_dx_X])+webRearSpar.qp_dx_X+ sum([webBottom.qp_dx_X])+webFrontSpar.qp_dx_X; %cell 1
421 Added: Fx = Fx + sum([webLowerNose.qp_dx_X])+ sum([webUpperNose.qp_dx_X]); %cell 2
422 Added: Fx
423 Added: Fz = sum([webTop.qp_dz_X])+webRearSpar.qp_dz_X+ sum([webBottom.qp_dz_X])+webFrontSpar.qp_dz_X; %cell 1
424 Added: Fz = Fz + sum([webLowerNose.qp_dz_X])+ sum([webUpperNose.qp_dz_X]); %cell 2
425 Added: Fz
426 Added:
427 Added: %check that q'*dz sums up to Vz
428 Added:
429 Added:
430 Added: Fx = sum([webTop.qp_dx_Z])+webRearSpar.qp_dx_Z+ sum([webBottom.qp_dx_Z])+webFrontSpar.qp_dx_Z; %cell 1
431 Added: Fx = Fx + sum([webLowerNose.qp_dx_Z])+ sum([webUpperNose.qp_dx_Z]); %cell 2
432 Added: Fx
433 Added: Fz = sum([webTop.qp_dz_Z])+webRearSpar.qp_dz_Z+ sum([webBottom.qp_dz_Z])+webFrontSpar.qp_dz_Z; %cell 1
434 Added: Fz = Fz + sum([webLowerNose.qp_dz_Z])+ sum([webUpperNose.qp_dz_Z]); %cell 2
435 Added: Fz
436 Added:
437 Added: %%
438 Added:
439 Added: % sum up the ds/t and q*ds/t to solve 2 equations, 2 unknowns
440 Added:
441 Added: % [A]*[q1s q2s] = B
442 Added:
443 Added: A11 = sum([webTop.dS_over_t])+webRearSpar.dS_over_t+ sum([webBottom.dS_over_t])+webFrontSpar.dS_over_t;
444 Added: A22 = sum([webLowerNose.dS_over_t])+ sum([webUpperNose.dS_over_t])+webFrontSparCell2.dS_over_t;
445 Added: A12 = -webFrontSpar.dS_over_t;
446 Added: A21 = -webFrontSparCell2.dS_over_t;
447 Added:
448 Added: B1_X = sum([webTop.q_dS_over_t_X])+webRearSpar.q_dS_over_t_X+ sum([webBottom.q_dS_over_t_X])+webFrontSpar.q_dS_over_t_X;
449 Added: B2_X = sum([webLowerNose.q_dS_over_t_X])+ sum([webUpperNose.q_dS_over_t_X])+webFrontSparCell2.q_dS_over_t_X;
450 Added: B1_Z = sum([webTop.q_dS_over_t_Z])+webRearSpar.q_dS_over_t_Z+ sum([webBottom.q_dS_over_t_Z])+webFrontSpar.q_dS_over_t_Z;
451 Added: B2_Z = sum([webLowerNose.q_dS_over_t_Z])+ sum([webUpperNose.q_dS_over_t_Z])+webFrontSparCell2.q_dS_over_t_Z;
452 Added:
453 Added: Amat = [A11 A12; A21 A22];
454 Added: Bmat_X = -[B1_X;B2_X];
455 Added: Bmat_Z = -[B1_Z;B2_Z];
456 Added:
457 Added: qs_X = inv(Amat)*Bmat_X;
458 Added: qs_Z = inv(Amat)*Bmat_Z;
459 Added:
460 Added:
461 Added:
462 Added: sum_2_a_q_X = sum([webTop.two_A_qprime_X])+webRearSpar.two_A_qprime_X+ sum([webBottom.two_A_qprime_X]); %cell 1 qprimes
463 Added: sum_2_a_q_X = sum_2_a_q_X + sum([webLowerNose.two_A_qprime_X])+ sum([webUpperNose.two_A_qprime_X]); %cell 2 qprimes
464 Added: sum_2_a_q_X = sum_2_a_q_X + 2*qs_X(1)*(sum([webTop.Area])+webRearSpar.Area+ sum([webBottom.Area]));
465 Added: sum_2_a_q_X = sum_2_a_q_X + 2*qs_X(2)*(sum([webLowerNose.Area])+ sum([webUpperNose.Area]));
466 Added:
467 Added: sum_2_a_q_Z = sum([webTop.two_A_qprime_Z])+webRearSpar.two_A_qprime_Z+ sum([webBottom.two_A_qprime_Z]); %cell 1 qprimes
468 Added: sum_2_a_q_Z = sum_2_a_q_Z + sum([webLowerNose.two_A_qprime_Z])+ sum([webUpperNose.two_A_qprime_Z]); %cell 2 qprimes
469 Added: sum_2_a_q_Z = sum_2_a_q_Z + 2*qs_Z(1)*(sum([webTop.Area])+webRearSpar.Area+ sum([webBottom.Area]));
470 Added: sum_2_a_q_Z = sum_2_a_q_Z + 2*qs_Z(2)*(sum([webLowerNose.Area])+ sum([webUpperNose.Area]));
471 Added:
472 Added: %shear center
473 Added: sc.posX = sum_2_a_q_Z / Vz + frontSpar*chord;
474 Added: sc.posZ = - sum_2_a_q_X / Vx;
475 Added:
476 Added:
477 Added: % now consider the torque representing shifting the load from the quarter
478 Added: % chord to the SC (need to check signs on these moments)
479 Added:
480 Added: torque_Z = Vz*(sc.posX - 0.25*chord);
481 Added: torque_X = -Vx*sc.posZ;
482 Added:
483 Added:
484 Added: Area1 = sum([webTop.Area]) + webRearSpar.Area + sum([webBottom.Area]);
485 Added: %check area
486 Added: Area1_check = get_int(frontSpar,backSpar,1)*chord^2 + get_int(frontSpar,backSpar,0)*chord^2;
487 Added:
488 Added: Area2 = sum([webLowerNose.Area]) + sum([webUpperNose.Area]);
489 Added: Area2_check = get_int(0,frontSpar,1)*chord^2 + get_int(0,frontSpar,0)*chord^2;
490 Added:
491 Added:
492 Added: %for twist equation (see excel spreadsheet example)
493 Added:
494 Added: q1t_over_q2t = (A22/Area2 + webFrontSpar.dS_over_t/Area1)/(A11/Area1 + webFrontSpar.dS_over_t/Area2);
495 Added:
496 Added: q2t = torque_X/(2*Area1*q1t_over_q2t + 2*Area2);
497 Added: q1t = q2t*q1t_over_q2t;
498 Added: qt_X = [q1t;q2t];
499 Added:
500 Added: q2t = torque_Z/(2*Area1*q1t_over_q2t + 2*Area2);
501 Added: q1t = q2t*q1t_over_q2t;
502 Added: qt_Z = [q1t;q2t];
503 Added:
504 Added:
505 Added:
506 Added: % --- - add up all shear flows: qtot = (qPrime + qs) + qt
507 Added:
508 Added:
509 Added:
510 Added:
511 Added: %--- insert force balance to check total shear flows ---
512 Added:
513 Added: % --- --
514 Added:
515 Added:
516 Added: %end
517 Added:
518 Added: sc
519 Added:
520 Added:
521 Added: %plotting airfoil cross-section
522 Added:
523 Added: xChord = 0:.01:1;
524 Added: xChord = xChord*chord;
525 Added: upperSurface = zeros(1,length(xChord));
526 Added: lowerSurface = zeros(1,length(xChord));
527 Added:
528 Added: for i=1:length(xChord)
529 Added: upperSurface(i) = get_z(xChord(i)/chord,1)*chord;
530 Added: lowerSurface(i) = get_z(xChord(i)/chord,0)*chord;
531 Added: end
532 Added:
533 Added: figure; hold on; axis equal; grid on;
534 Added: %plot(xChord,z_camber,'-')
535 Added: plot(xChord,upperSurface,'-k','linewidth',2)
536 Added: plot(xChord,lowerSurface,'-k','linewidth',2)
537 Added: plot([0 1],[0 0],'--k','linewidth',1)
538 Added:
539 Added:
540 Added: for i = 1:length(webTop)
541 Added: vecX = [frontSpar*chord webTop(i).xStart webTop(i).xEnd];
542 Added: vecZ = [0 webTop(i).zStart webTop(i).zEnd];
543 Added: fill(vecX,vecZ,[0.9 0.9 0.9])
544 Added: end
545 Added:
546 Added: for i = 1:length(webBottom)
547 Added: vecX = [frontSpar*chord webBottom(i).xStart webBottom(i).xEnd];
548 Added: vecZ = [0 webBottom(i).zStart webBottom(i).zEnd];
549 Added: fill(vecX,vecZ,[0.9 0.9 0.9])
550 Added: end
551 Added:
552 Added: for i = 1:length(webUpperNose)
553 Added: vecX = [frontSpar*chord webUpperNose(i).xStart webUpperNose(i).xEnd];
554 Added: vecZ = [0 webUpperNose(i).zStart webUpperNose(i).zEnd];
555 Added: fill(vecX,vecZ,[0.7 0.9 1.0])
556 Added: end
557 Added:
558 Added: for i = 1:length(webLowerNose)
559 Added: vecX = [frontSpar*chord webLowerNose(i).xStart webLowerNose(i).xEnd];
560 Added: vecZ = [0 webLowerNose(i).zStart webLowerNose(i).zEnd];
561 Added: fill(vecX,vecZ,[0.7 0.9 1.0])
562 Added: end
563 Added:
564 Added: vecX = [frontSpar*chord sparCaps(3).posX sparCaps(4).posX];
565 Added: vecZ = [0 sparCaps(3).posZ sparCaps(4).posZ];
566 Added: fill(vecX,vecZ,[0.9 0.9 0.9])
567 Added:
568 Added:
569 Added: sparCapSize = 18;
570 Added: stringerSize = 18;
571 Added: plot([sparCaps(1).posX sparCaps(2).posX],[sparCaps(1).posZ sparCaps(2).posZ],'-k','linewidth',2)
572 Added: plot([sparCaps(3).posX sparCaps(4).posX],[sparCaps(3).posZ sparCaps(4).posZ],'-k','linewidth',2)
573 Added: plot([sparCaps.posX],[sparCaps.posZ],'.b','markersize',sparCapSize)
574 Added: plot([topStringers.posX],[topStringers.posZ],'.r','markersize',stringerSize)
575 Added: plot([bottomStringers.posX],[bottomStringers.posZ],'.r','markersize',stringerSize)
576 Added: plot([noseTopStringers.posX],[noseTopStringers.posZ],'.r','markersize',stringerSize)
577 Added: plot([noseBottomStringers.posX],[noseBottomStringers.posZ],'.r','markersize',stringerSize)
578 Added: plot(centroid.posX,centroid.posZ,'.k','markerSize',18)
579 Added: plot(sc.posX,sc.posZ,'.g','markersize',18)