*ARCHIVED* development moved to aircraft-studio.
1
%wing shear flow
2
clear all;
3
close all;
4
5
Vx = 1; Vz = 1; My = 1; %test loads will be applied individually
6
7
8
%Ixz = -Ixz;
9
10
%define webs
11
12
%% web cell 1
13
14
%upper webs
15
numStringers = numTopStringers;
16
stringerGap = upperStringerGap;
17
webThickness = t_upper;
18
tempStringers = topStringers;
19
20
for i=1:(+1)
21
web().xStart = sparCaps(1).posX + stringerGap*(-1);
22
web().xEnd = sparCaps(1).posX + stringerGap*();
23
web().thickness = webThickness;
24
web().zStart = get_z(web()/,1)*chord;
25
web().zEnd = get_z(web()/,1)*chord;
26
if i==1
27
web().dp_area = sparCaps(1).area;
28
web().dP_X = 0;
29
web().dP_Z = 0;
30
web().qPrime_X = 0;
31
web().qPrime_Z = 0;
32
else
33
web().dp_area = tempStringers(-1).area;
34
dx = web().xStart-centroid.posX; dz = web().zStart-centroid.posZ;
35
web().dP_X = get_dp(,,,0,,,,web()); %just Vx
36
web().dP_Z = get_dp(,,0,,,,,web()); %just Vz
37
web().qPrime_X = web(-1).qPrime_X - web().dP_X;
38
web().qPrime_Z = web(-1).qPrime_Z - web().dP_Z;
39
end
40
tempInt = get_int(web()/,web()/,1)*chord^2; %integral of airfoil function
41
triangle1 = abs((web()-sparCaps(1))*web()/2);
42
triangle2 = abs((web()-sparCaps(1))*web()/2);
43
web().Area = tempInt + triangle1 - triangle2;
44
web().ds = get_ds(web()/,web()/,1)*chord;
45
web().dS_over_t = web().ds / web().thickness;
46
47
web().q_dS_over_t_X = web().qPrime_X * web().dS_over_t;
48
web().q_dS_over_t_Z = web().qPrime_Z * web().dS_over_t;
49
web().two_A_qprime_X = 2*web().Area*web().qPrime_X;
50
web().two_A_qprime_Z = 2*web().Area*web().qPrime_Z;
51
web().qp_dx_X = web().qPrime_X *(web()-web());
52
web().qp_dx_Z = web().qPrime_Z *(web()-web());
53
web().qp_dz_X = web().qPrime_X *(web()-web());
54
web().qp_dz_Z = web().qPrime_Z *(web()-web());
55
end
56
webTop = web;
57
web = [];
58
59
%rear spar
60
i=1;
61
web().xStart = sparCaps(3).posX;
62
web().xEnd = sparCaps(4).posX;
63
web().thickness = t_rearSpar;
64
web().zStart = sparCaps(3).posZ;
65
web().zEnd = sparCaps(4).posZ;
66
web().dp_area = sparCaps(3).area;
67
dx = web().xStart-centroid.posX; dz = web().zStart-centroid.posZ;
68
web().dP_X = get_dp(,,,0,,,,web());
69
web().dP_Z = get_dp(,,0,,,,,web());
70
web().qPrime_X = webTop(+1).qPrime_X - web().dP_X;
71
web().qPrime_Z = webTop(+1).qPrime_Z - web().dP_Z;
72
73
web().Area = (sparCaps(3)-sparCaps(1))*sparCaps(3).posZ/2 + ...
74
abs((sparCaps(3)-sparCaps(1))*sparCaps(4)/2);
75
web().ds = abs(sparCaps(3)-sparCaps(4));
76
web().dS_over_t = web().ds / web().thickness;
77
78
web().q_dS_over_t_X = web().qPrime_X * web().dS_over_t;
79
web().q_dS_over_t_Z = web().qPrime_Z * web().dS_over_t;
80
web().two_A_qprime_X = 2*web().Area*web().qPrime_X;
81
web().two_A_qprime_Z = 2*web().Area*web().qPrime_Z;
82
web().qp_dx_X = web().qPrime_X *(web()-web());
83
web().qp_dx_Z = web().qPrime_Z *(web()-web());
84
web().qp_dz_X = web().qPrime_X *(web()-web());
85
web().qp_dz_Z = web().qPrime_Z *(web()-web());
86
87
webRearSpar = web;
88
web = [];
89
90
91
%lower webs
92
numStringers = numBottomStringers;
93
stringerGap = lowerStringerGap;
94
webThickness = t_lower;
95
tempStringers = bottomStringers;
96
97
for i=1:(+1)
98
web().xStart = sparCaps(4).posX - stringerGap*(-1);
99
web().xEnd = sparCaps(4).posX - stringerGap*();
100
web().thickness = webThickness;
101
web().zStart = get_z(web()/,0)*chord;
102
web().zEnd = get_z(web()/,0)*chord;
103
dx = web().xStart-centroid.posX; dz = web().zStart-centroid.posZ;
104
if i==1
105
web().dp_area = sparCaps(4).area;
106
web().dP_X = get_dp(,,,0,,,,web());
107
web().dP_Z = get_dp(,,0,,,,,web());
108
web().qPrime_X = webRearSpar.qPrime_X - web().dP_X;
109
web().qPrime_Z = webRearSpar.qPrime_Z - web().dP_Z;
110
else
111
web().dp_area = tempStringers(-1).area;
112
web().dP_X = get_dp(,,,0,,,,web());
113
web().dP_Z = get_dp(,,0,,,,,web());
114
web().qPrime_X = web(-1).qPrime_X - web().dP_X;
115
web().qPrime_Z = web(-1).qPrime_Z - web().dP_Z;
116
end
117
118
tempInt = get_int(web()/,web()/,0)*chord^2; %integral of airfoil function
119
triangle2 = abs((web()-sparCaps(1))*web()/2);
120
triangle1 = abs((web()-sparCaps(1))*web()/2);
121
web().Area = tempInt + triangle1 - triangle2;
122
web().ds = get_ds(web()/,web()/,0)*chord;
123
web().dS_over_t = web().ds / web().thickness;
124
125
web().q_dS_over_t_X = web().qPrime_X * web().dS_over_t;
126
web().q_dS_over_t_Z = web().qPrime_Z * web().dS_over_t;
127
web().two_A_qprime_X = 2*web().Area*web().qPrime_X;
128
web().two_A_qprime_Z = 2*web().Area*web().qPrime_Z;
129
web().qp_dx_X = web().qPrime_X*(web()-web());
130
web().qp_dx_Z = web().qPrime_Z*(web()-web());
131
web().qp_dz_X = web().qPrime_X*(web()-web());
132
web().qp_dz_Z = web().qPrime_Z*(web()-web());
133
134
%web().radCurv = ... Example: get_curve(web(),web(),1)
135
end
136
webBottom = web;
137
web = [];
138
139
%front Spar
140
i=1;
141
web().xStart = sparCaps(2).posX;
142
web().xEnd = sparCaps(1).posX;
143
web().thickness = t_frontSpar;
144
web().zStart = sparCaps(2).posZ;
145
web().zEnd = sparCaps(1).posZ;
146
web().dp_area = sparCaps(2).area;
147
dx = web().xStart-centroid.posX; dz = web().zStart-centroid.posZ;
148
web().dP_X = get_dp(,,,0,,,,web());
149
web().dP_Z = get_dp(,,0,,,,,web());
150
web().qPrime_X = webBottom(+1).qPrime_X - web().dP_X;
151
web().qPrime_Z = webBottom(+1).qPrime_Z - web().dP_Z;
152
web().Area = 0;
153
web().ds = abs(sparCaps(2)-sparCaps(1));
154
web().dS_over_t = web().ds / web().thickness;
155
156
web().q_dS_over_t_X = web().qPrime_X * web().dS_over_t;
157
web().q_dS_over_t_Z = web().qPrime_Z * web().dS_over_t;
158
web().two_A_qprime_X = 2*web().Area*web().qPrime_X;
159
web().two_A_qprime_Z = 2*web().Area*web().qPrime_Z;
160
web().qp_dx_X = web().qPrime_X *(web()-web());
161
web().qp_dx_Z = web().qPrime_Z *(web()-web());
162
web().qp_dz_X = web().qPrime_X *(web()-web());
163
web().qp_dz_Z = web().qPrime_Z *(web()-web());
164
165
webFrontSpar = web;
166
web = [];
167
168
169
170
171
%% web cell 2
172
173
%lower nose webs
174
numStringers = numNoseBottomStringers;
175
stringerGap = lowerNoseStringerGap;
176
webThickness = t_lower_front;
177
tempStringers = noseBottomStringers;
178
179
for i=1:(+1)
180
web().xStart = sparCaps(2).posX - stringerGap*(-1);
181
web().xEnd = sparCaps(2).posX - stringerGap*();
182
web().thickness = webThickness;
183
web().zStart = get_z(web()/,0)*chord;
184
web().zEnd = get_z(web()/,0)*chord;
185
dx = web().xStart-centroid.posX; dz = web().zStart-centroid.posZ;
186
187
if i==1
188
web().dp_area = sparCaps(2).area;
189
web().dP_X = 0;
190
web().dP_Z = 0;
191
web().qPrime_X = 0;
192
web().qPrime_Z = 0;
193
else
194
web().dp_area = tempStringers(-1).area;
195
web().dP_X = get_dp(,,,0,,,,web());
196
web().dP_Z = get_dp(,,0,,,,,web());
197
web().qPrime_X = web(-1).qPrime_X - web().dP_X;
198
web().qPrime_Z = web(-1).qPrime_Z - web().dP_Z;
199
end
200
tempInt = get_int(web()/,web()/,0)*chord^2; %integral of airfoil function
201
triangle1 = abs((web()-sparCaps(2))*web()/2);
202
triangle2 = abs((web()-sparCaps(2))*web()/2);
203
web().Area = tempInt + triangle1 - triangle2;
204
web().ds = get_ds(web()/,web()/,0)*chord;
205
web().dS_over_t = web().ds / web().thickness;
206
207
web().q_dS_over_t_X = web().qPrime_X * web().dS_over_t;
208
web().q_dS_over_t_Z = web().qPrime_Z * web().dS_over_t;
209
web().two_A_qprime_X = 2*web().Area*web().qPrime_X;
210
web().two_A_qprime_Z = 2*web().Area*web().qPrime_Z;
211
web().qp_dx_X = web().qPrime_X *(web()-web());
212
web().qp_dx_Z = web().qPrime_Z *(web()-web());
213
web().qp_dz_X = web().qPrime_X *(web()-web());
214
web().qp_dz_Z = web().qPrime_Z *(web()-web());
215
216
%web().radCurv = ... Example: get_curve(web(),web(),1)
217
end
218
webLowerNose = web;
219
web = [];
220
221
%upper nose webs
222
numStringers = numNoseTopStringers;
223
stringerGap = upperNoseStringerGap;
224
webThickness = t_upper_front;
225
tempStringers = noseTopStringers;
226
227
for i=1:(+1)
228
web().xStart = stringerGap*(-1);
229
web().xEnd = stringerGap*();
230
web().thickness = webThickness;
231
web().zStart = get_z(web()/,1)*chord;
232
web().zEnd = get_z(web()/,1)*chord;
233
dx = web().xStart-centroid.posX; dz = web().zStart-centroid.posZ;
234
if i==1
235
web().dp_area = 0;
236
web().dP_X = 0;
237
web().dP_Z = 0;
238
web().qPrime_X = webLowerNose(+1).qPrime_X - web().dP_X;
239
web().qPrime_Z = webLowerNose(+1).qPrime_Z - web().dP_Z;
240
else
241
web().dp_area = tempStringers(-1).area;
242
web().dP_X = get_dp(,,,0,,,,web());
243
web().dP_Z = get_dp(,,0,,,,,web());
244
web().qPrime_X = web(-1).qPrime_X - web().dP_X;
245
web().qPrime_Z = web(-1).qPrime_Z - web().dP_Z;
246
end
247
tempInt = get_int(web()/,web()/,1)*chord^2; %integral of airfoil function
248
triangle2 = abs((web()-sparCaps(2))*web()/2);
249
triangle1 = abs((web()-sparCaps(2))*web()/2);
250
web().Area = tempInt + triangle1 - triangle2;
251
web().ds = get_ds(web()/,web()/,1)*chord;
252
web().dS_over_t = web().ds / web().thickness;
253
254
web().q_dS_over_t_X = web().qPrime_X * web().dS_over_t;
255
web().q_dS_over_t_Z = web().qPrime_Z * web().dS_over_t;
256
web().two_A_qprime_X = 2*web().Area*web().qPrime_X;
257
web().two_A_qprime_Z = 2*web().Area*web().qPrime_Z;
258
web().qp_dx_X = web().qPrime_X *(web()-web());
259
web().qp_dx_Z = web().qPrime_Z *(web()-web());
260
web().qp_dz_X = web().qPrime_X *(web()-web());
261
web().qp_dz_Z = web().qPrime_Z *(web()-web());
262
263
end
264
webUpperNose = web;
265
web = [];
266
267
268
%front Spar
269
i=1;
270
web().xStart = sparCaps(1).posX;
271
web().xEnd = sparCaps(2).posX;
272
web().thickness = t_frontSpar;
273
web().zStart = sparCaps(1).posZ;
274
web().zEnd = sparCaps(2).posZ;
275
web().dp_area = sparCaps(1).area;
276
dx = web().xStart-centroid.posX; dz = web().zStart-centroid.posZ;
277
278
web().dP_X = get_dp(,,,0,,,,web());
279
web().dP_Z = get_dp(,,0,,,,,web());
280
web().qPrime_X = webUpperNose(+1).qPrime_X - web().dP_X;
281
web().qPrime_Z = webUpperNose(+1).qPrime_Z - web().dP_Z;
282
web().Area = 0;
283
web().ds = abs(sparCaps(1)-sparCaps(2));
284
web().dS_over_t = web().ds / web().thickness;
285
web().q_dS_over_t_X = web().qPrime_X * web().dS_over_t;
286
web().q_dS_over_t_Z = web().qPrime_Z * web().dS_over_t;
287
web().two_A_qprime_X = 2*web().Area*web().qPrime_X;
288
web().two_A_qprime_Z = 2*web().Area*web().qPrime_Z;
289
web().qp_dx_X = web().qPrime_X *(web()-web());
290
web().qp_dx_Z = web().qPrime_Z *(web()-web());
291
web().qp_dz_X = web().qPrime_X *(web()-web());
292
web().qp_dz_Z = web().qPrime_Z *(web()-web());
293
294
webFrontSparCell2 = web;
295
web = [];
296
297
298
%check that q'*dx sums up to Vx
299
300
Fx = sum([webTop.qp_dx_X])+webRearSpar.qp_dx_X+ sum([webBottom.qp_dx_X])+webFrontSpar.qp_dx_X; %cell 1
301
Fx = Fx + sum([webLowerNose.qp_dx_X])+ sum([webUpperNose.qp_dx_X]); %cell 2
302
Fx
303
Fz = sum([webTop.qp_dz_X])+webRearSpar.qp_dz_X+ sum([webBottom.qp_dz_X])+webFrontSpar.qp_dz_X; %cell 1
304
Fz = Fz + sum([webLowerNose.qp_dz_X])+ sum([webUpperNose.qp_dz_X]); %cell 2
305
Fz
306
307
%check that q'*dz sums up to Vz
308
309
310
Fx = sum([])+webRearSpar.qp_dx_Z+ sum([])+webFrontSpar.qp_dx_Z; %cell 1
311
Fx = Fx + sum([])+ sum([]); %cell 2
312
Fx
313
Fz = sum([])+webRearSpar.qp_dz_Z+ sum([])+webFrontSpar.qp_dz_Z; %cell 1
314
Fz = Fz + sum([])+ sum([]); %cell 2
315
Fz
316
317
%%
318
319
% sum up the ds/t and q*ds/t to solve 2 equations, 2 unknowns
320
321
% []*[q2s] = B
322
323
A11 = sum([dS_over_t])+webRearSpar.dS_over_t+ sum([dS_over_t])+webFrontSpar.dS_over_t;
324
A22 = sum([dS_over_t])+ sum([dS_over_t])+webFrontSparCell2.dS_over_t;
325
A12 = -webFrontSpar.dS_over_t;
326
A21 = -webFrontSparCell2.dS_over_t;
327
328
B1_X = sum([])+webRearSpar.q_dS_over_t_X+ sum([])+webFrontSpar.q_dS_over_t_X;
329
B2_X = sum([])+ sum([])+webFrontSparCell2.q_dS_over_t_X;
330
B1_Z = sum([])+webRearSpar.q_dS_over_t_Z+ sum([])+webFrontSpar.q_dS_over_t_Z;
331
B2_Z = sum([])+ sum([])+webFrontSparCell2.q_dS_over_t_Z;
332
333
Amat = [;A22];
334
Bmat_X = -[;];
335
Bmat_Z = -[;];
336
337
qs_X = inv()*Bmat_X;
338
qs_Z = inv()*Bmat_Z;
339
340
341
342
sum_2_a_q_X = sum([])+webRearSpar.two_A_qprime_X+ sum([]); %cell 1 qprimes
343
sum_2_a_q_X = sum_2_a_q_X + sum([])+ sum([]); %cell 2 qprimes
344
sum_2_a_q_X = sum_2_a_q_X + 2*qs_X(1)*(sum([])++sum([]));
345
sum_2_a_q_X = sum_2_a_q_X + 2*qs_X(2)*(sum([])+sum([]));
346
347
sum_2_a_q_Z = sum([])+webRearSpar.two_A_qprime_Z+ sum([]); %cell 1 qprimes
348
sum_2_a_q_Z = sum_2_a_q_Z + sum([])+ sum([]); %cell 2 qprimes
349
sum_2_a_q_Z = sum_2_a_q_Z + 2*qs_Z(1)*(sum([])++sum([]));
350
sum_2_a_q_Z = sum_2_a_q_Z + 2*qs_Z(2)*(sum([])+sum([]));
351
352
%shear center
353
sc.posX = sum_2_a_q_Z / Vz + frontSpar*chord;
354
sc.posZ = - sum_2_a_q_X / Vx;
355
356
357
% now consider the torque representing shifting the load from the quarter
358
% chord to the SC(moments)
359
360
torque_Z = Vz*(-0.25*);
361
torque_X = -Vx*sc.posZ;
362
363
364
Area1 = sum([]) + webRearSpar.Area + sum([]);
365
%check area
366
Area1_check = get_int(,,1)*chord^2 + get_int(,,0)*chord^2;
367
368
Area2 = sum([]) + sum([]);
369
Area2_check = get_int(0,,1)*chord^2 + get_int(0,,0)*chord^2;
370
371
372
%for twist equation(example)
373
374
q1t_over_q2t = (/+dS_over_t/)/(/+dS_over_t/);
375
376
q2t = torque_X/(2**+2*);
377
q1t = q2t*q1t_over_q2t;
378
qt_X = [;];
379
380
q2t = torque_Z/(2**+2*);
381
q1t = q2t*q1t_over_q2t;
382
qt_Z = [;];
383
384
385
386
% --- - add up all shear flows: qtot = (+) + qt
387
388
389
390
391
%--- insert force balance to check total shear flows ---
392
393
% --- --
394
395
396
%end
397
398
sc
399
400
401
%plotting airfoil cross-section
402
403
xChord = 0:.01:1;
404
xChord = xChord*chord;
405
upperSurface = zeros(1,length());
406
lowerSurface = zeros(1,length());
407
408
for i=1:length()
409
upperSurface() = get_z(xChord()/,1)*chord;
410
lowerSurface() = get_z(xChord()/,0)*chord;
411
end
412
413
figure; hold on; axis equal; grid on;
414
%plot(,,'-')
415
plot(,,'-k','linewidth',2)
416
plot(,,'-k','linewidth',2)
417
plot([01],[00],'--k','linewidth',1)
418
419
420
for i = 1:length()
421
vecX = [*webTopwebTop];
422
vecZ = [0webTopwebTop];
423
fill(,,[0.90.90.9])
424
end
425
426
for i = 1:length()
427
vecX = [*webBottomwebBottom];
428
vecZ = [0webBottomwebBottom];
429
fill(,,[0.90.90.9])
430
end
431
432
for i = 1:length()
433
vecX = [*webUpperNosewebUpperNose];
434
vecZ = [0webUpperNosewebUpperNose];
435
fill(,,[0.70.91.0])
436
end
437
438
for i = 1:length()
439
vecX = [*webLowerNosewebLowerNose];
440
vecZ = [0webLowerNosewebLowerNose];
441
fill(,,[0.70.91.0])
442
end
443
444
vecX = [*sparCaps3sparCaps4];
445
vecZ = [0sparCaps3sparCaps4];
446
fill(,,[0.90.90.9])
447
448
449
sparCapSize = 18;
450
stringerSize = 18;
451
plot([sparCaps1sparCaps2],[sparCaps1sparCaps2],'-k','linewidth',2)
452
plot([sparCaps3sparCaps4],[sparCaps3sparCaps4],'-k','linewidth',2)
453
plot([],[],'.b','markersize',)
454
plot([],[],'.r','markersize',)
455
plot([],[],'.r','markersize',)
456
plot([],[],'.r','markersize',)
457
plot([],[],'.r','markersize',)
458
plot(,,'.k','markerSize',18)
459
plot(,,'.g','markersize',18)
460