User:Eas4200c.f08.aeris.krammer/HW7

Minor edit. Formatting of picture sizing edited. Content remains the same.
Eas4200c.f08.aeris.mcrae 23:30, 12 December 2008 (UTC)

Unsymmetric Thin-Walled Cross Sections

Figure 1: Asymmetric thin-walled cross section with stringers of varying areas.

Figure 2: Asymmetric thin-walled contour cross section.

We know that for a symmetric thin-walled cross section that $\ q_{s}=-\int \int _{A_{s}}^{}{{\frac {d\sigma _{xx}}{dx}}dA}$

For unsymmetric thin-walled cross sections under bidirectional bending, the bending stress $\left(\sigma _{xx}\right)$ must be taken into account. Doing so, we arrive at the following equation for sheer flow.

$\ q_{s}=-\left(k_{y}V_{y}-k_{yz}V_{z}\right)Q_{z}+\left(k_{z}V_{z}-k_{yz}V_{y}\right)Q_{y}$

We can begin our analysis of the figures to the right using the Mean Value Theorem [1] to find the first moments of area around the z and y axes.

$\ Q_{z}=\int \int _{A_{s}}^{}{ydA}={\bar {y}}A$

$\ Q_{y}=\int \int _{A_{s}}^{}{zdA}={\bar {z}}A$

For the case of stringers of different areas, as shown in Figure 1, the area to be used in this equation is $\ A=\sum _{i=1}^{4}{A_{i}}$. For the case shown in Figure 2, the area used in this equation is $\ A=A_{\hat {s}}+A_{{\hat {s}},L}$.

Using this, as well as equations from the previous homework, we can break down analysis of unsymmetric thin-walled cross sections, such as the two shown on the right in Figure 1 and Figure 2, into the following four steps.

1. Find $\left({\bar {y}},{\bar {z}}\right)$.

2. Find $I_{y}$, $I_{z}$, and $I_{yz}$.

3. Find $k_{y}$, $k_{z}$, and $k_{yz}$.

4. Follow path $s$ to find $q_{12}$, $q_{23}$, and $q_{34}$.

Shear Flow in a Single and Multi-Cell Section

There are different scenarios that can be analyzed to obtain a better understanding of shear flow in sections. The scenarios are divided into Single-Cell sections and Multi-Cell sections. In a Single-Cell section, two more scenarios can be analyzed. One without stringers and another with stringers. This can also be done for Multi-Cell sections.

The following plan was used to abbreviate and guide the analysis of shear flow through each different scenario:

S) Single-Cell sections

S1) Single-Cell section without stringers (most basic analysis)
S2) Single-Cell section with stringers

M) Multi-Cell sections

M1) Multi-Cell section without stringers
M2) Multi-Cell section with stringers

To analyze S1, the cell is looked at as a thin-walled cross-section. The shear flow through the section must be found.

Shear Flow Across S1 (HW student contribution)
$R^{z}=R_{BA}^{z}+R_{AB}^{z}=-qA'B'=-R_{BA}^{z}\Rightarrow R^{z}=0$

To analyze the shear flow in S2, the problem is divided into two parts using the superposition principal. The first part is to obtain the shear flow in a cell without stringers. The second part is to analyze a cell with stringers and applying a cut in one of its shear paths. The shear flows from each part is then added together to obtain the shear flow in S2.

The shear flow within each section, between stringers, are not equal to each other but are constant.

Therefore shear flow of S2 is also constant.

To conclude the analysis, solve for the shear flow in the second part and take the moment about a point selected on the plane of the section.

Multi-Cell Sections

Using superposition, we are able to use the following equation when solving for shear flow:

                         qij = q + qijN

where,q_ij, q, and q_ij^N are true shear flow, closed cell constant shear/flow, and open cell piecewise constant shear flow, respectively. This property can be applied to a multi-cell section with or without stringers.

Case 1: Without Stringers

R^z = ΣR_i^z where R^z is the resultant of the z-axis. This equation explains a very useful property for this case: R_i^z = 0 so, R^z = 0 This resultant can be expanded to any number of cells.

Case 2: With Stringers

Using the equation from the previous case, we can apply it to this case by

                         Rz1 = ΣRiz1 = 0

This is true for this case because when the summation is expanded, we get: R^z = R^z₁ + R^z₂. Each term is the resultant of P, resultant of P₁ and resultant of P₂ respectively. For this case, the resultants of P and P₂ are equal to each other which causes P₁ = 0.

To reinforce the use of this property, let us do an example and calculate the equilibrium for each stringer in P₂.

STRINGER CALCULATIONS (EULER-CUT PRINCIPLE)

ΣFx = 0 = $\int \limits _{A3}$[σxx(x + dx) - σxx(x)]dA3 [-${\tilde {q}}23$ - ${\tilde {q}}43$ + ${\tilde {q}}31$]dx This series, when solved for ${\tilde {q}}31$ is ${\tilde {q}}31={\tilde {q}}23+{\tilde {q}}43+q^{3}$ where $q^{3}$ is the contribution to shear flow by stringer 3. Solving, we get q3 = $\int \limits _{A3}$dσxxdxdA3. Summing the additional forces of transverse shear (Vy and Vz) we get a final response of q3 = -(kyVy - kyzVz)Qz3 - (kzVz - kyzVy)Qy3 where, Qz3 = y3A3 (Y-coordinate of the centroid and Area of stringer 3) and Qy3 = z3A3. ${\tilde {q}}23={\tilde {q}}12-{\tilde {q}}24+q^{2}$. Since ${\tilde {q}}12$ and ${\tilde {q}}24$ take place at a cut, they are equivalent to 0 making ${\tilde {q}}23=q^{2}$. Applying the same principles as the previous stringer we get: q2 = -(kyVy - kyzVz)Qz2 - (kzVz - kyzVy)Qy2 where, Qz2 = y2A2 and Qy2 = z2A2. ${\tilde {q}}43={\tilde {q}}24-{\tilde {q}}41+q^{4}$. Since ${\tilde {q}}24$ and ${\tilde {q}}41$ take place at a cut, they are also equivalent to 0 making ${\tilde {q}}43=q^{4}$. Applying the principle once again, we get: q4 = -(kyVy - kyzVz)Qz4 - (kzVz - kyzVy)Qy4 where, Qz4 = y4A4 and Qy4 = z4A4. Now that we know what all the terms are for our equation of ${\tilde {q}}31={\tilde {q}}23+{\tilde {q}}43+q^{3}$, we can plug them in and use superposition to give us: $q_{23}={\tilde {q}}23+q_{1}$ $q_{31}={\tilde {q}}31+q_{1}$ $q_{24}={\tilde {q}}24+q_{2}$ $q_{43}={\tilde {q}}43+q_{2}$ $q_{41}={\tilde {q}}41+q_{3}$ where q1, q2, and q3 are unknown. Since we have 3 unknowns, we need 3 equations to solve. We can use: 1) Moment Equation - Pick the most convenient point that will make the calculations easy (lines of action) 2) θ1 = θ2 (Compatibility Equation) (i.e. Rate of Twist) 3) θ2 = θ3 (Compatibility Equation) (i.e. Rate of Twist) =

Equilibrium in Stringers

Equilibrium needs to be obtained in all the stingers when calculating shear flow. There are 2 ways to compute the shear flow: one method has been shown earlier, and the second way to compute the shear flow is the consequence of this first method.

Stringer with shear flow in and out.

We can analyze what happens of we cut the walls and isolate one stringer.

Shear flow in a cell can be analyzed by making a cut so it can be treated as an open section with an existing constant shear flow. We can deduce the following relation from the figure with an isolated stringer due to the cuts:

Cell with cut walls - isolated stringer.

Isolated stringer.

$\displaystyle {\tilde {q}}_{23}={\tilde {q}}_{31}={\tilde {q}}_{34}=0$

We can also see from the equilibrium equation that

$\displaystyle {\tilde {q}}_{23}={\tilde {q}}_{31}={\tilde {q}}_{34}+q^{(3)}$

From the above equations we can see that $\displaystyle q^{(3)}$ must be 0, which cannot be true.
Therefore we can say that in making cuts no part of the section should be isolated.

Solving P2

The resulting equation for shear flow would be

$\displaystyle 0=q^{(1)}+q^{(2)}+q^{(4)}$

However, this relationship is not true, and it is not even possible. The following explains why.

$\displaystyle 0=q^{(1)}+q^{(2)}+q^{(3)}+q^{(4)}$

Or this relationship can be rewritten as

$\displaystyle 0=\sum _{4}^{e=1}{}q^{(e)}$

Therefore,

$\displaystyle q^{(1)}+q^{(2)}+q^{(4)}=-q^{(3)}\neq 0$

We can prove that the sum of all shear flows is zero

$\displaystyle q^{(e)}=n_{z}Q_{z}^{(e)}+n_{y}Q_{y}^{(e)}$

Where,

$\displaystyle n_{z}:=-(k_{y}V_{y}-k_{yz}V_{z})$

$\displaystyle n_{y}:=-(k_{z}V_{z}-k_{yz}V_{y})$

Defining the first area moment of inertia

Therefore the sum of all shear flows is

$\displaystyle \sum _{4}^{e=1}{}q^{(e)}=n_{z}\sum _{4}^{e=1}{}Q_{z}^{(e)}+n_{y}\sum _{4}^{e=1}{}Q_{y}^{(e)}$

However, $\displaystyle \sum _{4}^{e=1}{}Q_{y}^{(e)}=0$ and $\displaystyle \sum _{4}^{e=1}{}Q_{z}^{(e)}=0$

Since Q is the first area moment of inertia, computed respect to the coordinate system which has its origin centered at y-z.

Therefore,

$\displaystyle Q({\hat {z}})=\int _{A({\hat {z}})}^{}{}zdA$

And

$\displaystyle Q({\hat {y}})=0=z_{c}\int _{A}^{}{}dA$

The relationship above is due to the mean value theorem.

Shear Flow

In order to plot the buckling shape, we express the coefficients $\left\{{\begin{matrix}C_{22}&C_{13}&C_{31}&C_{33}\end{matrix}}\right\}$ in terms of $C_{11}$. Using an aspect ratio $\left(\vartheta \right)$ equal to 1.5 we can break our analysis down to the following four steps.

1. Find $\lambda$.

$\displaystyle \lambda =\left[{\frac {\vartheta ^{4}}{81(1+\vartheta ^{2})^{4}}}\left\{1+{\frac {81}{625}}+{\frac {81}{25}}\left({\frac {1+\vartheta ^{2}}{1+9\vartheta ^{2}}}\right)^{2}+{\frac {81}{25}}\left({\frac {1+\vartheta ^{2}}{9+\vartheta ^{2}}}\right)^{2}\right\}\right]^{1/2}=0.0288$

2. Evaluate $K_{5x5}$ numerically.

Failed to parse (syntax error): {\displaystyle \ \left[ \begin{array}{lllll} \frac{\lambda (1 + \vartheta^2)^2}{\vartheta^2} & \frac{4}{9} & 0 & 0 & 0 \\ \frac{4}{9} & \frac{16 \lambda (1 + \vartheta^2)^2}{\vartheta^2} & - \frac{4}{5} & - \frac{4}{5} & \frac{36}{25} \\ 0 & - \frac{4}{5} & \frac{\lambda (1 + 9 \vartheta^2)^2}{\vartheta^2} & 0 & 0 \\ 0 & - \frac{4}{5} & 0 & \frac{\lambda (9 + \vartheta^2)^2}{\vartheta^2} & 0 \\ 0 & \frac{36}{25} & 0 & 0 & \frac{\lambda (9 + 9 \vartheta^2)^2}{\vartheta^2} \end{array} \right] \left\{ \begin{array}{l} C_{11} \\ C_{22} \\ C_{13} \\ C_{31} \\ C_{33} \end{array} \right\} = \left\{ \begin{array}{l} 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{array} \right\} }

3. Truncate $K_{5x5}$ into $K_{4x4}$.

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "http://localhost:6011/en.wikiversity.org/v1/":): {\displaystyle \ \mathbf \bar K _{4 \times 4} = \left[ \begin{array}{llll} \frac{16 \lambda (1 + \vartheta^2)^2}{\vartheta^2} & - \frac{4}{5} & - \frac{4}{5} & \frac{36}{25} \\ - \frac{4}{5} & \frac{\lambda (1 + 9 \vartheta^2)^2}{\vartheta^2} & 0 & 0 \\ - \frac{4}{5} & 0 & \frac{\lambda (9 + \vartheta^2)^2}{\vartheta^2} & 0 \\ \frac{36}{25} & 0 & 0 & \frac{\lambda (9 + 9 \vartheta^2)^2}{\vartheta^2} \end{array} \right] \left\{ \begin{array}{l} C_{22} \\ C_{13} \\ C_{31} \\ C_{33} \end{array} \right\} = \left\{ \begin{array}{l} \frac{-4 C_{11}}{9} \\ 0 \\ 0 \\ 0 \end{array} \right\} }

4. Solve for $\left\{{\begin{matrix}C_{22}&C_{13}&C_{31}&C_{33}\end{matrix}}\right\}$ in terms of $C_{11}$, noting that $K^{-1}$ is the inverse of $K$.

${\begin{Bmatrix}C_{22}\\C_{13}\\C_{31}\\C_{33}\end{Bmatrix}}=K^{-1}{\begin{Bmatrix}-{\frac {4}{9}}C_{11}\\0\\0\\0\end{Bmatrix}}$

$u_{z}=C_{11}sin({\frac {\pi x}{a}})*sin({\frac {\pi y}{b}})+C_{22}sin({\frac {2\pi x}{a}})*sin({\frac {2\pi y}{b}})+C_{31}sin({\frac {\pi x}{a}})*sin({\frac {3\pi y}{b}})+C_{31}sin({\frac {3\pi x}{a}})*sin({\frac {\pi y}{b}})+C_{33}sin({\frac {3\pi x}{a}})*sin({\frac {3\pi y}{b}})$

Bidirectional Bending Continued

Let us continue to examine the beam with an arbitrarily-shaped cross section.
As shown in the figure, we assume that the x-axis coincides with the centroidal axis. We also assume that no torsion is present since the loads (P) pass through the shear center.

Bidirectional bending of a bar with an arbitrarily shaped cross-section

Under this bidirectional bending the displacement is a function of x, y and z variables. The displacement equations can be defined as:

$\displaystyle u=u_{0}(x)+z\psi _{y}(x)+y\psi _{y}(x)$

$\displaystyle v=v_{0}(x)$

$\displaystyle w=w_{0}(x)$

Where $\displaystyle u_{0}$ denotes the displacement in the x direction, $\displaystyle v_{0}$ denotes the displacement in the y direction and $\displaystyle w_{0}$ denotes the displacement in the z direction. The subscript 0 refers to the neutral axis. $\displaystyle \psi _{y}$ and $\displaystyle \psi _{z}$ are rotations about the cross section y and z axes respectively.

The corresponding strains are

$\varepsilon _{xx}={\frac {\partial u}{\partial x}}={\frac {du_{0}}{dx}}+z{\frac {d\psi _{y}}{dx}}+y{\frac {d\psi _{z}}{dx}}$

$\gamma _{xy}={\frac {\partial v}{\partial x}}+{\frac {\partial u}{\partial y}}={\frac {dv_{0}}{dx}}+\psi _{z}$

$\gamma _{xz}={\frac {\partial w}{\partial x}}+{\frac {\partial u}{\partial z}}={\frac {dw_{0}}{dx}}+\psi _{y}$

We can assume that $\displaystyle \gamma _{xy}=\gamma _{xz}=0$ . This yields to
$\displaystyle \gamma _{z}=-{\frac {dv_{0}}{dx}}$
$\displaystyle \gamma _{y}=-{\frac {dw_{0}}{dx}}$

Substituting these results back to the bending strain equation we get

$\varepsilon _{xx}={\frac {du_{0}}{dx}}-y{\frac {d^{2}v_{0}}{dx^{2}}}-z{\frac {d^{2}w_{0}}{dx^{2}}}$

Since, $\displaystyle {\frac {du_{0}}{dx}}=0$ and $\displaystyle N_{x}=0$, we can simplify the equation above:

$\varepsilon _{xx}=-y{\frac {d^{2}v_{0}}{dx^{2}}}-z{\frac {d^{2}w_{0}}{dx^{2}}}$

Therefore, we can define the bending moments about the y and z axes

$M_{y}=\int _{A}^{}{}\int z\sigma _{xx}dA=-E\int _{A}^{}{}\int (yz{\frac {dv_{0}}{dx^{2}}}+z^{2}{\frac {dw_{0}}{dx^{2}}})dA=-EI_{yz}{\frac {d^{2}v_{0}}{dx^{2}}}-EI_{y}{\frac {d^{2}w_{0}}{dx^{2}}}$

$M_{z}=\int _{A}^{}{}\int y\sigma _{xx}dA==-EI_{z}{\frac {d^{2}v_{0}}{dx^{2}}}-EI_{yz}{\frac {d^{2}w_{0}}{dx^{2}}}$

Where the inertia moments are defined such
The moment of inertia about the y-axis
$\displaystyle I_{y}=\int _{A}^{}{}\int z^{2}dA$
The moment of inertia about the z-axis
$\displaystyle I_{z}=\int _{A}^{}{}\int y^{2}dA$
And the product of inertia is
$\displaystyle I_{yz}=\int _{A}^{}{}\int yzdA$

The moment equations above can be expressed in matrix form using the following notations for the curvatures: $\displaystyle X_{y}={d^{2}v_{0}}{dx^{2}}$ and $\displaystyle X_{z}={d^{2}w_{0}}{dx^{2}}$

${\begin{Bmatrix}M_{y}\\M_{z}\end{Bmatrix}}={\begin{vmatrix}I_{y}&I_{yz}\\I_{yz}&I_{z}\end{vmatrix}}{\begin{Bmatrix}-EX_{z}\\-EX_{y}\end{Bmatrix}}$

Also, strain can be expressed in a similar manner using the curvature relations above.

$\varepsilon _{xx}=-yX_{y_{}}-zX_{z}={\begin{bmatrix}z&y\end{bmatrix}}{\begin{Bmatrix}-X_{z}\\-X_{y}\end{Bmatrix}}$

$\sigma _{xx}=E\varepsilon _{xx}=E{\begin{bmatrix}z&y\end{bmatrix}}{\begin{Bmatrix}-X_{z}\\-X_{y}\end{Bmatrix}}=E{\begin{bmatrix}z&y\end{bmatrix}}{\frac {1}{E}}I^{-1}{\begin{Bmatrix}M_{y}\\M_{z}\end{Bmatrix}}$

The inverse of the 2 by 2 matrix is

$I^{-1}={\frac {1}{D}}{\begin{vmatrix}I_{z}&-I_{yz}\\-I_{yz}&I_{y}\end{vmatrix}}$

And the determinant D is defined as $\displaystyle D=I_{y}I_{z}-(I_{yz})^{2}$

We can directly derive the martix form of the shear flow.
We know that the general expression for shear flow is

$q=-\int _{A}^{}{}{\frac {d\sigma _{xx}}{dx}}dA$

We also need a relationship between the sear and the moment

${\begin{matrix}d\sigma _{xx}\\dx\end{matrix}}={\begin{bmatrix}z&y\end{bmatrix}}I^{-1}{\begin{Bmatrix}{\frac {dM_{y}}{dx}}\\{\frac {dM_{z}}{dx}}\end{Bmatrix}}$

From here we can conclude that

$V_{y}={\frac {dM_{y}}{dx}},V_{z}={\frac {dM_{z}}{dx}}$
and

$Q_{y}=\int _{A}^{}{}zdA,Q_{z}=\int _{A}^{}{}ydA$

Matlab HW7: Determining the True Shear Flow in Different-Celled NACA Airfoils

     The problem statement in this homework was to determine the true shear flows of two different cases of the NACA airfoil subject to shear forces. Part 1 was to determine these shear flows in a single-cell airfoil and part 2 was to determine these shear flows in a double-cell airfoil. The airfoils in both cases have four stringers: 2 positioned at a quarter of the chord length and 2 positioned at three quarters of the chord length. In each case, the problem was solved using superposition of the shear flows in a case analyzed without stringers and a case with stringers. The equation of moment and the compatibility equation were also used to find the true shear flows. The code in which was developed to solve this problem is shown in the collapsible boxes below.

     For part 1, the code starts out with finding the values ky, kz, and kyz by using the moment of inertias obtained by previous homework assignments. Then, it goes about solving the P2 problem as posed in lecture. The trailing edge is the portion of the skin that is designated to be cut. Here, Qy and Qz are calculated at each stringer node and the equilibrium equations of the shear flows at each stringer are written. This solves for all the q tilda's in P2. After these shear flows are obtained, it is necessary to find a the shear flow of the two cells (q_11 and q_21) solving for the P1 problem as posed in lecture. The indices in the notation are described as follows: the first subscript is the cell number, and the second subscript is the part of the homework the shear flow corresponds to. To do this, the moment equation is used and the moment was taken about point D. To find the moments contributed by the shear flows found in P2, the code was written to sweep from stringer to stringer that had a non-zero shear flow. At each sweep, an infinitesimal length of the skin was calculated and multiplied by the shear flow through that skin portion to obtain a force. Once that force was determined, the perpendicular distance from point D to the line of action of the force was obtained and then multiplied to the force to get an infinitesimally small moment. All the moments were added together to get the moment produced by each shear flow from P2. Along with the compatibility equation, q_11 and q_21 are obtained. Back to superposition, the shear flows of P1 and P2 are added together to get the true shear flows.

     For part 2, the code uses the same moments as part 1 because nothing has changed. However, the change is in the moment equation. The moment equation now has three unknowns (q_12 and q_22 and q_32) so it alone is not sufficient in solving for the P1 shear flows. This is where the compatibility equation comes in. We know that the twist angle of the first cell is equal to the twist angle of the second and third cells. The moment and compatibility equations are arranged to that it has coefficients multiplied to q_12 and q_22 and q_32. Then, the equations were arranged in matrix form and solved for in Matlab. Then, superposition was used to find the true shear flow in the three-cell case.

Matlab Code to Determine True Shear Flows (HW student contribution)

function airfoil clc clear all close all %{ m = input('1st digit of the NACA airfoil: '); p = input('2nd digit of the NACA airfoil: '); t_a = input('3rd and 4th digits of the NACA airfoil: '); c = input('Chord length in meters: '); ns = input('Number of segments: '); %} % Location of reference point P_0y = 0; P_0z = 0; %Applied shear forces V_y=10000; %N V_z=5000; %N % Applied bending moments M_y = - 1250; % N.m M_z = + 500; % N.m % Stringer areas A_B = 2/10000; % m^2 A_E = 2/10000; % m^2 A_F = 1/10000; % m^2 A_H = 1/10000; % m^2 % Stringer thinckness t_str = .005; % m % Skin thickness and web thickness t_s = .002; % m t_w = .003; % m % Debug values c = .5; ns = 180; m = 2; p = 4; t_a =15; % Determine coordinates of the airfoil contour [y_u,y_l,z_u,z_l]=naca_plot(c,m,p,t_a,ns); % Find the index locations of the upper and lower contour % for 1/4 c and 3/4 c [garb,iB] = min(abs(y_u-c/4)); [garb,iE] = min(abs(y_l-c/4)); [garb,iF] = min(abs(y_u-3*c/4)); [garb,iH] = min(abs(y_l-3*c/4)); % Assign leading edge (L) and trailing edge (T) indicies iL = 1; iT = ns; % Define points B = [c/4; z_u(iB)]; D = [c/4; 0]; E = [c/4; z_l(iE)]; F = [3*c/4; z_u(iF)]; G = [3*c/4; 0]; H = [3*c/4; z_l(iH)]; % ============= OLD CODE ================== % Calculate area of cell 1 % Abar_1 = A_upper + A_lower Abar_1 = calc_area(D,y_u,z_u,iB,-1,iL) + calc_area(D,y_l,z_l,iL,1,iE); % Calculate area of cell 2 % This is done in two parts % 1 - At point E, sweep forward from F to B % 2 - At point F, sweep backwards from E to H Abar_2 = calc_area(E,y_u,z_u,iF,-1,iB) + calc_area(F,y_l,z_l,iE,1,iH); % Calculate area of cell 3 % Abar_3 = A_upper + A_lower Abar_3 = calc_area(G,y_u,z_u,iT,-1,iF) + calc_area(G,y_l,z_l,iT,-1,iH); Abar = Abar_1+Abar_2+Abar_3; % Calculate the length of each contour segment len_LB = calc_length(y_u,z_u,iL,iB); len_BF = calc_length(y_u,z_u,iB,iF); len_FT = calc_length(y_u,z_u,iF,iT); len_LE = calc_length(y_l,z_l,iL,iE); len_EH = calc_length(y_l,z_l,iE,iH); len_HT = calc_length(y_l,z_l,iH,iT); len_BE = z_u(iB) - z_l(iE); len_FH = z_u(iF) - z_l(iH); perim = len_LB + len_BF + len_FT + len_LE + len_EH + len_HT; % Calculate centroid of the skin cnt = 0; for i=1:ns-1 cnt = cnt + [(y_u(i)+y_u(i+1))* sqrt((z_u(i+1)-z_u(i))^2+(y_u(i+1)-y_u(i))^2)+(y_l(i)+y_l(i+1))* sqrt((z_l(i+1)-z_l(i))^2+(y_l(i+1)-y_l(i))^2); ... (z_u(i)+z_u(i+1))*sqrt((z_u(i+1)-z_u(i))^2+(y_u(i+1)-y_u(i))^2)+(z_l(i)+z_l(i+1))*sqrt((z_l(i+1)-z_l(i))^2+(y_l(i+1)-y_l(i))^2)]; end skin_cent = cnt / (2*(perim)); % Calculate the centroid of the spar webs spar_cent = (len_BE*(B+E) + len_FH*(F+H))/(2*(len_BE + len_FH)); % Single-cell J_singlecell = 4*(Abar)^2/(len_LB+len_BF+len_FT+len_LE+len_EH+len_HT)*t_s; % Multi-cell c1 = ((len_LB+len_LE)/t_s + len_BE/t_w)/Abar_1; c2 = len_BE/(t_w*Abar_1); c3 = ((len_BF+len_EH)/t_s + len_BE/t_w + len_FH/t_w)/Abar_2; c4 = len_BE/(t_w*Abar_2); c5 = len_FH/(t_w*Abar_2); c6 = ((len_FT+len_HT)/t_s + len_FH/t_w)/Abar_3; c7 = len_FH/(t_w*Abar_3); C1 = (c2*c5+c2*c6-c2*c5+c6*c3)/(c1*c5+c1*c6+c4*c6); C2 = (c1*c3+c1*c7-c4*c2+c4*c7)/(c1*c5+c1*c6+c4*c6); J_multicell = 4*(Abar_1*C1 + Abar_2 + Abar_3*C2)/(c3 - c4*C1 - c5*C2); %============= CASE 1: Only stringers =========== % Calculate stringer centroid A_tot = A_B + A_E + A_F + A_H; str_cent = (B*A_B + E*A_E + H*A_F + F*A_H)/A_tot; % Calculate the second area moment of inertia I_22 = (A_B*(B(2)-str_cent(2))^2) + (A_E*(E(2)-str_cent(2))^2) + (A_H*(H(2)-str_cent(2))^2) + (A_F*(F(2)-str_cent(2))^2); I_33 = (A_B*(B(1)-str_cent(1))^2) + (A_E*(E(1)-str_cent(1))^2) + (A_H*(H(1)-str_cent(1))^2) + (A_F*(F(1)-str_cent(1))^2); I_23 = (A_B*(B(2)-str_cent(2))*(B(1)-str_cent(1))) + (A_E*(E(2)-str_cent(2))*(E(1)-str_cent(1))) + (A_H*(H(2)-str_cent(2))*(H(1)-str_cent(1))) + (A_F*(F(2)-str_cent(2))*(F(1)-str_cent(1))); % Calculate the angle of the neutral axis alpha = atan((I_22*M_z - I_23*M_y)/(I_33*M_y-I_23*M_z)); % Slope of neutral axis (-alpha) beta = - alpha; % Generate two points to plot the neutral axis y_na = [-.2*c, 1.2*c]; z_na = beta * (y_na - str_cent(1)) + str_cent(2); % Calculate the bending stress at each point sig_11(1) = (I_22*M_z-I_23*M_y)/(I_22*I_33-I_23^2)*B(1) + (I_33*M_y-I_23*M_z)/(I_22*I_33-I_23^2)*B(2); sig_11(2) = (I_22*M_z-I_23*M_y)/(I_22*I_33-I_23^2)*E(1) + (I_33*M_y-I_23*M_z)/(I_22*I_33-I_23^2)*E(2); sig_11(3) = (I_22*M_z-I_23*M_y)/(I_22*I_33-I_23^2)*F(1) + (I_33*M_y-I_23*M_z)/(I_22*I_33-I_23^2)*F(2); sig_11(4) = (I_22*M_z-I_23*M_y)/(I_22*I_33-I_23^2)*H(1) + (I_33*M_y-I_23*M_z)/(I_22*I_33-I_23^2)*H(2); [garb,ind] = max(abs(sig_11)); switch ind case 1 sig_max = (I_22*M_z-I_23*M_y)/(I_22*I_33-I_23^2)*B(1) + (I_33*M_y-I_23*M_z)/(I_22*I_33-I_23^2)*B(2); sig_max_y = B(1); sig_max_z = B(2); case 2 sig_max = (I_22*M_z-I_23*M_y)/(I_22*I_33-I_23^2)*E(1) + (I_33*M_y-I_23*M_z)/(I_22*I_33-I_23^2)*E(2); sig_max_y = E(1); sig_max_z = E(2); case 3 sig_max = (I_22*M_z-I_23*M_y)/(I_22*I_33-I_23^2)*F(1) + (I_33*M_y-I_23*M_z)/(I_22*I_33-I_23^2)*F(2); sig_max_y = F(1); sig_max_z = F(2); case 4 sig_max = (I_22*M_z-I_23*M_y)/(I_22*I_33-I_23^2)*H(1) + (I_33*M_y-I_23*M_z)/(I_22*I_33-I_23^2)*H(2); sig_max_y = H(1); sig_max_z = H(2); end % Calculate the bending moments with ultimate yield stress of steel sig_u = 1860e6; M_ratio = M_y/M_z; M_uy = sig_u/((I_22/M_ratio - I_23)/(I_22*I_33 - I_23^2)*sig_max_y + (I_33 - I_23/M_ratio)/(I_22*I_33 - I_23^2)*sig_max_z); M_uz = sig_u/((I_22 - I_23*M_ratio)/(I_22*I_33 - I_23^2)*sig_max_y + (I_33*M_ratio - I_23)/(I_22*I_33 - I_23^2)*sig_max_z); % ============= CASE 2: Complete analysis ============= % Calculate the centroid of the entire airfoil cent = (perim*t_s*skin_cent + (len_BE + len_FH)*t_w*spar_cent + A_tot*str_cent)/(perim*t_s + (len_BE + len_FH)*t_w + A_tot); % Add the contribution of the skin to the moments of inertia (I += A*d^2) for i=1:ns-1 I_22 = I_22 + sqrt((z_u(i+1)-z_u(i))^2+(y_u(i+1)-y_u(i))^2) *t_s * ((z_u(i)+z_u(i+1))/2)^2 + ... sqrt((z_l(i+1)-z_l(i))^2+(y_l(i+1)-y_l(i))^2) * t_s * ((z_l(i)+z_l(i+1))/2)^2; I_33 = I_33 + sqrt((z_u(i+1)-z_u(i))^2+(y_u(i+1)-y_u(i))^2) *t_s * ((y_u(i)+y_u(i+1))/2)^2 + ... sqrt((z_l(i+1)-z_l(i))^2+(y_l(i+1)-y_l(i))^2) * t_s * ((y_l(i)+y_l(i+1))/2)^2; I_23 = I_23 + sqrt((z_u(i+1)-z_u(i))^2+(y_u(i+1)-y_u(i))^2) *t_s * (z_u(i)+z_u(i+1))*(y_u(i)+y_u(i+1))/4 + ... sqrt((z_l(i+1)-z_l(i))^2+(y_l(i+1)-y_l(i))^2) * t_s * (z_l(i)+z_l(i+1))*(y_l(i)+y_l(i+1))/4; end % Add the contribution of the spar webs to the moments of inertia I_22 = I_22 + t_w*(len_BE+len_FH)*spar_cent(2)^2; I_33 = I_33 + t_w*(len_BE+len_FH)*spar_cent(1)^2; I_23 = I_23 + t_w*(len_BE+len_FH)*spar_cent(1)*spar_cent(2); % Calculate the angle of the neutral axis alpha = atan((I_22*M_z - I_23*M_y)/(I_33*M_y-I_23*M_z)); % Slope of neutral axis (-alpha) beta = - alpha; % Generate two points to plot the neutral axis y_na = [-.2*c, 1.2*c]; z_na = beta * (y_na - cent(1)) + cent(2); % Calculate the bending stress at each point sig_11(1) = (I_22*M_z-I_23*M_y)/(I_22*I_33-I_23^2)*B(1) + (I_33*M_y-I_23*M_z)/(I_22*I_33-I_23^2)*B(2); sig_11(2) = (I_22*M_z-I_23*M_y)/(I_22*I_33-I_23^2)*E(1) + (I_33*M_y-I_23*M_z)/(I_22*I_33-I_23^2)*E(2); sig_11(3) = (I_22*M_z-I_23*M_y)/(I_22*I_33-I_23^2)*F(1) + (I_33*M_y-I_23*M_z)/(I_22*I_33-I_23^2)*F(2); sig_11(4) = (I_22*M_z-I_23*M_y)/(I_22*I_33-I_23^2)*H(1) + (I_33*M_y-I_23*M_z)/(I_22*I_33-I_23^2)*H(2); [garb,ind] = max(abs(sig_11)); switch ind case 1 sig_max = (I_22*M_z-I_23*M_y)/(I_22*I_33-I_23^2)*B(1) + (I_33*M_y-I_23*M_z)/(I_22*I_33-I_23^2)*B(2); sig_max_y = B(1); sig_max_z = B(2); case 2 sig_max = (I_22*M_z-I_23*M_y)/(I_22*I_33-I_23^2)*E(1) + (I_33*M_y-I_23*M_z)/(I_22*I_33-I_23^2)*E(2); sig_max_y = E(1); sig_max_z = E(2); case 3 sig_max = (I_22*M_z-I_23*M_y)/(I_22*I_33-I_23^2)*F(1) + (I_33*M_y-I_23*M_z)/(I_22*I_33-I_23^2)*F(2); sig_max_y = F(1); sig_max_z = F(2); case 4 sig_max = (I_22*M_z-I_23*M_y)/(I_22*I_33-I_23^2)*H(1) + (I_33*M_y-I_23*M_z)/(I_22*I_33-I_23^2)*H(2); sig_max_y = H(1); sig_max_z = H(2); end % Calculate the bending moments with ultimate yield stress of steel sig_u = 1860e6; M_ratio = M_y/M_z; M_uy = sig_u/((I_22/M_ratio - I_23)/(I_22*I_33 - I_23^2)*sig_max_y + (I_33 - I_23/M_ratio)/(I_22*I_33 - I_23^2)*sig_max_z); M_uz = sig_u/((I_22 - I_23*M_ratio)/(I_22*I_33 - I_23^2)*sig_max_y + (I_33*M_ratio - I_23)/(I_22*I_33 - I_23^2)*sig_max_z); % ================= end HW5 ====================== temp =0; for i=iB:iF sig_BF(i) = (I_22*M_z-I_23*M_y)/(I_22*I_33-I_23^2)*y_u(i) + (I_33*M_y-I_23*M_z)/(I_22*I_33-I_23^2)*z_u(i); temp = temp+sig_BF(i); end sig_BF_av = temp/(iF-iB+1); temp=0; for i=iE:iH sig_EH(i) = (I_22*M_z-I_23*M_y)/(I_22*I_33-I_23^2)*y_l(i) + (I_33*M_y-I_23*M_z)/(I_22*I_33-I_23^2)*z_l(i); temp = temp+sig_EH(i); end sig_EH_av = temp/(iH-iE+1); % figure(2) % hold on % subplot(2,1,1),plot(y_u(iB:iF),sig_BF(iB:iF),'k-',y_u(iB:iF),sig_BF_av,'b:') % %plot (y_u(iB:iF),sig_BF_av,'b:') % legend('sigma_x_x from B-F','sigma_a_v_e') % title('Plot of sig_x_x vs. distance in top and bottom panels') % hold off % xlabel('Distance (m)') % ylabel('sigma_x_x (Pa)') % % hold on % subplot(2,1,2), plot(y_u(iE:iH),sig_EH(iE:iH),'k-',y_u(iE:iH),sig_EH_av,'b:') % legend('sigma_x_x from E-H','sigma_a_v_e') % hold off % xlabel('Distance (m)') % ylabel('sigma_x_x (Pa)') % ================= Plot the results ================ % z-values for spar webs used to plot spars spar1_z = [z_u(iB), z_l(iE)]; spar2_z = [z_u(iF), z_l(iH)]; % Plot the airfoil contour and the spar webs % figure(1) % hold on % plot(y_u,z_u,'k-',y_l,z_l,'k-') % plot([c/4, c/4],spar1_z,'b:') % plot([3*c/4, 3*c/4],spar2_z,'b:') % plot(str_cent(1),str_cent(2),'r+') % plot(cent(1),cent(2),'*') % plot(B(1),B(2),'ko',E(1),E(2),'ko',H(1),H(2),'ko',F(1),F(2),'ko','MarkerFaceColor','k') % plot(y_na,z_na,'-.') % plot(sig_max_y,sig_max_z, 'r*') % hold off % title(['Profile of NACA 2415 Airfoil (c = ' num2str(c) ', ns = ' num2str(ns) ')']) % text(B(1)-.018*c,B(2)+.03*c, 'B') % text(E(1)-.018,E(2)-.03*c, 'E') % text(F(1)-.018,F(2)+.03*c, 'F') % text(H(1)-.018,H(2)-.03*c, 'H') % text(.2*c,.35*c,['Maximum bending stress = ' num2str(sig_max) 'Pa']) % text(.2*c,.30*c,'Moments of inertia components (I_2_2,I_3_3,I_2_3):') % text(.2*c,.26*c,[num2str(I_22) ' m^4 ' num2str(I_33) ' m^4 ' num2str(I_23) ' m^4 ']) % axis([-.1*c 1.1*c -.4*c .4*c]) % ylabel('z') % xlabel('y') %========================== HW7 part1 ========================== %Double cell NACA Airfoil. Assuming no spar webs BE and FH. k_y=I_22/(I_22*I_33-I_23^2); k_z=I_33/(I_22*I_33-I_23^2); k_yz=I_23/(I_22*I_33-I_23^2); % Solving P2 (Qy and Qz are calculated using the centroid of the entire % airfoil.) %analyze stringer F QyFB=(F(2)-cent(2))*A_F; QzFB=(F(1)-cent(1))*A_F; q_tildFB=-(k_y*V_y-k_yz*V_z)*QzFB-(k_z*V_z-k_yz*V_y)*QyFB; %N/m %analyze stringer B QyBE=(B(2)-cent(2))*A_B; QzBE=(B(1)-cent(1))*A_B; q_tildBE=-(k_y*V_y-k_yz*V_z)*QzBE-(k_z*V_z-k_yz*V_y)*QyBE+q_tildFB; %N/m q_tildEB=0; %N/m %analyze stringer E QyEH=(E(2)-cent(2))*A_E; QzEH=(E(1)-cent(1))*A_E; q_tildEH=-(k_y*V_y-k_yz*V_z)*QzEH-(k_z*V_z-k_yz*V_y)*QyEH+q_tildBE-q_tildEB; %N/m %analyze stringer H q_tildHF=0; %Taking the moment about point D (which is on the line of action of V_y and %V_z momentFB=0; for i=iF:-1:iB %calculating a very small moment for each partitioned length and corresponding shear flow from F to B len=sqrt((y_u(i)- y_u(i+1))^2 + (z_u(i)-z_u(i+1))^2); force=q_tildFB*len; m=(z_u(i+1)-z_u(i))/(y_u(i+1)-y_u(i)); slope=-1/m; y_int=(slope-m)^-1*(z_u(i)-m*y_u(i)-D(2)+slope*D(1)); z_int=slope*(y_int-D(1))+D(2); momentarm=sqrt((y_int-D(1))^2+(z_int-D(2))^2); dmoment=force*momentarm; momentFB=momentFB+dmoment; end momentFB; momentBL=0; for i=iB:-1:1 %calculating a very small moment for each partitioned length and corresponding shear flow from B to L len=sqrt((y_u(i)- y_u(i+1))^2 + (z_u(i)-z_u(i+1))^2); force=q_tildBE*len; m=(z_u(i+1)-z_u(i))/(y_u(i+1)-y_u(i)); slope=-1/m; y_int=(slope-m)^-1*(z_u(i)-m*y_u(i)-D(2)+slope*D(1)); z_int=slope*(y_int-D(1))+D(2); momentarm=sqrt((y_int-D(1))^2+(z_int-D(2))^2); dmoment=force*momentarm; momentBL=momentBL+dmoment; end momentBL; momentLE=0; for i=1:iL %calculating a very small moment for each partitioned length and corresponding shear flow from L to E len=sqrt((y_l(i)- y_l(i+1))^2 + (z_l(i)-z_l(i+1))^2); force=q_tildBE*len; m=(z_l(i+1)-z_l(i))/(y_l(i+1)-y_l(i)); slope=-1/m; y_int=(slope-m)^-1*(z_l(i)-m*y_l(i)-D(2)+slope*D(1)); z_int=slope*(y_int-D(1))+D(2); momentarm=sqrt((y_int-D(1))^2+(z_int-D(2))^2); dmoment=force*momentarm; momentLE=momentLE+dmoment; end momentLE; momentEH=0; for i=iE:iH %calculating a very small moment for each partitioned length and corresponding shear flow from E to H len=sqrt((y_l(i)- y_l(i+1))^2 + (z_l(i)-z_l(i+1))^2); force=q_tildBE*len; m=(z_l(i+1)-z_l(i))/(y_l(i+1)-y_l(i)); slope=-1/m; y_int=(slope-m)^-1*(z_l(i)-m*y_l(i)-D(2)+slope*D(1)); z_int=slope*(y_int-D(1))+D(2); momentarm=sqrt((y_int-D(1))^2+(z_int-D(2))^2); dmoment=force*momentarm; momentEH=momentEH+dmoment; end momentEH; %Calculating using the moment equation and compatibility equation. Solving %the system of linear equations by matrix algebra. Syscoef=[2*Abar_1,2*Abar_2;(-len_BE/(t_w*Abar_2)-(len_LB+len_LE)/(t_s*Abar_1)-len_BE/(t_w*Abar_1)),(len_BE/(t_w*Abar_2)+len_EH/(t_s*Abar_2)+len_FH/(t_w*Abar_2)+len_BF/(t_s*Abar_2)+len_BE/(t_w*Abar_1))]; Knowns=[-(momentFB+momentBL+momentLE+momentEH);0]; Shearflows=Syscoef^-1*Knowns; Knowns=[-(momentFB+momentBL+momentLE+momentEH);0]; Shearflows=Syscoef^-1*Knowns; q_11=Shearflows(1); %N/m, shear flow in cell 1 for part 2 q_21=Shearflows(2); %N/m, shear flow in cell 2 for part 2 %Back to superposition. disp('These are the true shear flows for the airfoil for part 1. Answers are in N/m.') disp('A negative sign denotes that the shear flow is in the opposite direction of the assumed direction') q_FB1=q_21+q_tildFB q_BE1=q_11+q_tildBE q_EB1=q_11-q_21+q_tildEB q_EH1=q_21+q_tildEH q_HF1=q_21 %========================== HW7 part2 ========================== %3 cell NACA Airfoil. Assuming spar webs BE and FH. % Solving P2 (Qy and Qz are calculated using the centroid of the entire % airfoil.) We will solve a system of 3 equations to find q_12, q_22, and % q_32 specifcally using the moment equation and the compatibility equation. %analyze stringer F QyFB=(F(2)-cent(2))*A_F; QzFB=(F(1)-cent(1))*A_F; q_tildFB=-(k_y*V_y-k_yz*V_z)*QzFB-(k_z*V_z-k_yz*V_y)*QyFB; %N/m %analyze stringer B QyBE=(B(2)-cent(2))*A_B; QzBE=(B(1)-cent(1))*A_B; q_tildBE=-(k_y*V_y-k_yz*V_z)*QzBE-(k_z*V_z-k_yz*V_y)*QyBE+q_tildFB; %N/m q_tildEB=0; %N/m %analyze stringer E QyEH=(E(2)-cent(2))*A_E; QzEH=(E(1)-cent(1))*A_E; q_tildEH=-(k_y*V_y-k_yz*V_z)*QzEH-(k_z*V_z-k_yz*V_y)*QyEH+q_tildBE; %N/m %analyze stringer H q_tildHF=0; %Taking the moment about point D (which is on the line of action of V_y and %V_z momentFB=0; for i=iF:-1:iB %calculating a very small moment for each partitioned length and corresponding shear flow from F to B len=sqrt((y_u(i)- y_u(i+1))^2 + (z_u(i)-z_u(i+1))^2); force=q_tildFB*len; m=(z_u(i+1)-z_u(i))/(y_u(i+1)-y_u(i)); slope=-1/m; y_int=(slope-m)^-1*(z_u(i)-m*y_u(i)-D(2)+slope*D(1)); z_int=slope*(y_int-D(1))+D(2); momentarm=sqrt((y_int-D(1))^2+(z_int-D(2))^2); dmoment=force*momentarm; momentFB=momentFB+dmoment; end momentFB; momentBL=0; for i=iB:-1:1 %calculating a very small moment for each partitioned length and corresponding shear flow from B to L len=sqrt((y_u(i)- y_u(i+1))^2 + (z_u(i)-z_u(i+1))^2); force=q_tildBE*len; m=(z_u(i+1)-z_u(i))/(y_u(i+1)-y_u(i)); slope=-1/m; y_int=(slope-m)^-1*(z_u(i)-m*y_u(i)-D(2)+slope*D(1)); z_int=slope*(y_int-D(1))+D(2); momentarm=sqrt((y_int-D(1))^2+(z_int-D(2))^2); dmoment=force*momentarm; momentBL=momentBL+dmoment; end momentBL; momentLE=0; for i=1:iL %calculating a very small moment for each partitioned length and corresponding shear flow from L to E len=sqrt((y_l(i)- y_l(i+1))^2 + (z_l(i)-z_l(i+1))^2); force=q_tildBE*len; m=(z_l(i+1)-z_l(i))/(y_l(i+1)-y_l(i)); slope=-1/m; y_int=(slope-m)^-1*(z_l(i)-m*y_l(i)-D(2)+slope*D(1)); z_int=slope*(y_int-D(1))+D(2); momentarm=sqrt((y_int-D(1))^2+(z_int-D(2))^2); dmoment=force*momentarm; momentLE=momentLE+dmoment; end momentLE; momentEH=0; for i=iE:iH %calculating a very small moment for each partitioned length and corresponding shear flow from E to H len=sqrt((y_l(i)- y_l(i+1))^2 + (z_l(i)-z_l(i+1))^2); force=q_tildBE*len; m=(z_l(i+1)-z_l(i))/(y_l(i+1)-y_l(i)); slope=-1/m; y_int=(slope-m)^-1*(z_l(i)-m*y_l(i)-D(2)+slope*D(1)); z_int=slope*(y_int-D(1))+D(2); momentarm=sqrt((y_int-D(1))^2+(z_int-D(2))^2); dmoment=force*momentarm; momentEH=momentEH+dmoment; end momentEH; %Calculating using the moment equation and compatibility equation. Solving %the system of linear equations by matrix algebra. Syscoef=[2*Abar_1,2*Abar_2,2*Abar_3; (-len_BE/(t_w*Abar_2)-(len_LB+len_LE)/(t_s*Abar_1)-len_BE/(t_w*Abar_1)),... (len_BE/(t_w*Abar_2)+len_EH/(t_s*Abar_2)+len_FH/(t_w*Abar_2)+len_BF/(t_s*Abar_2)+len_BE/(t_w*Abar_1)),... (-len_FH/(t_w*Abar_2));... (-len_BE/(t_w*Abar_2)),(len_BF/(t_s*Abar_2)+len_BE/(t_w*Abar_2)+len_EH/(t_s*Abar_2)+len_FH/(t_w*Abar_2)+len_FH/(t_w*Abar_3)),... (-len_FH/(t_w*Abar_2)+(len_HT+len_FT)/(t_s*Abar_3)+len_FH/(t_w*Abar_3));]; Knowns=[-(momentFB+momentBL+momentLE+momentEH);0;0]; Shearflows=Syscoef^-1*Knowns; q_12=Shearflows(1); %N/m, shear flow in cell 1 for part 2 q_22=Shearflows(2); %N/m, shear flow in cell 2 for part 2 q_32=Shearflows(3); %N/m, shear flow in cell 3 for part 2 %Back to superposition. disp('These are the true shear flows for the airfoil for part 2. Answers are in N/m.') disp('A negative sign denotes that the shear flow is in the opposite direction of the assumed direction') q_FB2=q_22+q_tildFB q_BE2=q_11+q_tildBE q_EB2=q_11-q_22 q_EH2=q_22+q_tildEH q_HF2=-q_22+q_32+q_tildHF end %{ naca_plot: This function will return the upper and lower coordinates of the airfoil %} function [y_u,y_l,z_u,z_l]=naca_plot(c,m,p,t_a,ns) p = p/10; % position of maximum camber in tenths of the chord length c m = m/100; % maximum camber in percentage of chord length c t_a = t_a/100; % maximum airfoil thickness in percentage of the chord c % Create an array of ns elements between 0 and 1 y = linspace(0,1,ns); % This loop creates the mean camber line contour (z) and slope (dz) for i = 1: ns % ns+1 because we start at 0 and go to ns if y(i) <= p z_c(i) = (m/p^2) * (2*p*y(i) - y(i)^2); dz_c(i) = (m/p^2) * (2*p - 2*y(i)); else z_c(i) = (m/(1 - p)^2) * ((1 - 2*p)+ 2*p*y(i) - y(i)^2); dz_c(i) = (m/(1-p)^2) * (2*p - 2*y(i)); end end % calculate the distance from the camber line and the tangential angle %z_t = abs(t_a/0.2 * (.2969*sqrt(y) - .1260*y. - .3516*(y.)^2 + .2843*(y.)^3 - .1015*(y.)^4)); z_t = abs(5.*t_a.*(.2969.*sqrt(y)-.1260.*y-.3516.*y.^2 + .2843.*y.^3 - .1015*y.^4)); theta = atan(dz_c); % Calculate the upper and lower coordinates of the foil and scale to chord % length (c) y_u = (y - z_t.*sin(theta))*c; z_u = (z_c + z_t.*cos(theta))*c; y_l = (y + z_t.*sin(theta))*c; z_l = (z_c - z_t.*cos(theta))*c; end %{ area: This function calculates the area of shape given the contour coordinates Inputs - A : Reference point y : vector of y-coordinates z : vector of z-coordinates i_start : starting index i_end : stopping index step : increment for each loop (usually 1 or -1) %} function area = calc_area(A,y,z,i_start,step,i_end) area = 0; for i = i_start:step:i_end-step B = [y(i); z(i)]; C = [y(i+step); z(i+step)]; area = area + .5*cross([B-A;0],[C-A;0]); % Must be 3x3 vector end area = norm(real(area)); end %{ perimeter: This function calculates the length of a contour given the coordinates Inputs - y : vector of y-coordinates z : vector of z-coordinates i_start : starting index i_end : stopping index %} function perimeter = calc_length(y,z,i_start,i_end) perimeter = 0; for i = i_start:i_end-1 len = sqrt((y(i)- y(i+1))^2 + (z(i)-z(i+1))^2); perimeter = perimeter + len; end end

Matlab Code Output(HW student contribution)

These are the true shear flows for the airfoil for part 1. Answers are in N/m. A negative sign denotes that the shear flow is in the opposite direction of the assumed direction q_FB1 = 5.1989e+003 q_BE1 = -1.1903e+004 q_EB1 = -878.8807 q_EH1 = 4.4195e+003 q_HF1 = 9.4520e+003 These are the true shear flows for the airfoil for part 2. Answers are in N/m. A negative sign denotes that the shear flow is in the opposite direction of the assumed direction q_FB2 = 5.6237e+003 q_BE2 = -1.1903e+004 q_EB2 = -1.3037e+003 q_EH2 = 4.8444e+003 q_HF2 = -1.3959e+004

Next we analyzed buckling from shear flow. The critical buckling stress is analyzed two different ways. First we approximate the critical buckling by truncating the equation below to two variables, $C_{11}$ and $C_{22}$.

$\displaystyle \left[{\begin{array}{lllll}{\frac {\lambda (1+\vartheta ^{2})^{2}}{\vartheta ^{2}}}&{\frac {4}{9}}&0&0&0\\{\frac {4}{9}}&{\frac {16\lambda (1+\vartheta ^{2})^{2}}{\vartheta ^{2}}}&-{\frac {4}{5}}&-{\frac {4}{5}}&{\frac {36}{25}}\\0&-{\frac {4}{5}}&{\frac {\lambda (1+9\vartheta ^{2})^{2}}{\vartheta ^{2}}}&0&0\\0&-{\frac {4}{5}}&0&{\frac {\lambda (9+\vartheta ^{2})^{2}}{\vartheta ^{2}}}&0\\0&{\frac {36}{25}}&0&0&{\frac {\lambda (9+9\vartheta ^{2})^{2}}{\vartheta ^{2}}}\end{array}}\right]\left\{{\begin{array}{l}C_{11}\\C_{22}\\C_{13}\\C_{31}\\C_{33}\end{array}}\right\}=\left\{{\begin{array}{l}0\\0\\0\\0\\0\end{array}}\right\}$

becomes,

$\displaystyle \left[{\begin{array}{ll}{\frac {\lambda (1+\vartheta ^{2})^{2}}{\vartheta ^{2}}}&{\frac {4}{9}}\\{\frac {4}{9}}&{\frac {16\lambda (1+\vartheta ^{2})^{2}}{\vartheta ^{2}}}\end{array}}\right]\left\{{\begin{array}{l}C_{11}\\C_{22}\end{array}}\right\}=\left\{{\begin{array}{l}0\\0\end{array}}\right\}$

Since we are truncating the original equation, we lose some accuracy. By keeping the original 5x5 matrix and solving for the eigenvalues, we obtain the second equation for critical buckling stress. Therefore, the two equations for critical buckling stress become:

Failed to parse (syntax error): {\displaystyle \displaystyle (\sigma_{xy})_{cr}^{\star,2} = \pm k_c^{\star,2} \frac{\pi^2 D}{b^2 h} \ , \ k_c^{\star,2} = \frac{\pi^2}{32} \frac{9 (1 + \vartheta^2)^2}{\vartheta^3} }

And

Failed to parse (syntax error): {\displaystyle \displaystyle (\sigma_{xy})_{cr}^{\star,5} = \pm k_c^{\star,5} \frac{\pi^2 D}{b^2 h} \ , \ k_c^{\star,5} = \frac{\pi^2}{32} \frac{1}{\lambda \vartheta} }

By plotting the two equations with the parameters from the airfoil panel FB, we can see the error that occurs by truncating. It also becomes evident that the error increases at larger values of $a/b$.

If we add the plot for a clamped boundary condition we see that the critical buckling stress increases as expected. From the plot it is evident that a value of $K=4.1$ is not possible; hence, this value was an error in the experiment. A larger plot is shown below.

Matlab Code to analyze buckling

% a) a=sym('[lambda*(1+v^2)^2/v^2 4/9; 4/9 16*lambda*(1+v^2)^2/v^2]'); ad=det(a); solve(ad,'lambda') % b) ns = 100; E = 72e9; % Pa nu = .33; b = 0.2520; % m h = .002; % m v = linspace(.5,2,ns); m = 1; D = E*h^3/(12*(1-nu^2)); for i=1:ns sig_c2(i) = (pi^2*9*(1+v(i)^2)^2/(32*v(i)^3)) *pi^2*D/(b^2*h); end % c) for i=1:ns lambda(i) = (v(i)^4/(81*(1+v(i)^2)^4) * (1+81/625+81/25*((1+v(i)^2)/(1+9*v(i)^2))^2+81/25*((1+v(i)^2)/(9+v(i)^2))^2))^.5; sig_c5(i) = (pi^2/(32*lambda(i)*v(i))) *pi^2*D/(b^2*h); end figure(1) plot (v,sig_c2,v,sig_c5); title('Plot of sigma_c_r^2 and sigma_c_r^5 vs. a/b for panel BF') xlabel('a/b') ylabel('sigma_c_r') legend('sigma_c_r^2','sigma_c_r^5','location','best') % d) K = [12.7 9.5 7.38 4.1]; sig_c = K*E/(1-nu^2)*(h/b)^2; figure(2) hold on plot (v,sig_c2,v,sig_c5); %plot (v,sig_c5,'r'); plot ([1,2,999], sig_c(1:3),'r-'); plot([0,999],[sig_c(4),sig_c(4)],'k--') hold off axis([0.5 4.5 0 16e7]) title('Plot of sigma_c_r^2 and sigma_c_r^5 and clamped sigma_c_r vs. a/b for panel BF') xlabel('a/b') ylabel('sigma_c_r') legend('sigma_c_r^2','sigma_c_r^5','sigma_c_r clamped','sigma_c_r K=4.1','location','best') figure(3) hold on plot ([1,2,99], sig_c(1:3),'r-'); plot([0,99],[sig_c(4),sig_c(4)],'k--') hold off axis([0.5 20 0 16e7]) title('Plot of sigma_c_r (clamped condition) vs. a/b for panel BF') xlabel('a/b') ylabel('sigma_c_r') legend('sigma_c_r','sigma_c_r K=4.1','location','best')

Relating the work above to the NACA airfoil, we can calculate the shear stress through the panel BF with the calculated shear flow.

$N_{xy}^{}=\sigma _{xy}h$

$\sigma _{xy}=N_{xy}^{}/h$

Since the thickness, $h$, of the panel is given as 2 mm, the shear stress is calculated as follows:

$\sigma _{xy}={\frac {5623.7(N/m)}{0.002(m)}}=2,811,850(Pa)$

For large values of a/b, this value approaches the critical buckling predicted above. Since this is an airfoil, it is safe to assume $a/b$ will be much smaller than one. For this assumption, the stresses on the airfoil will not induce buckling of the panels.

E-learning versus MediaWiki

     From previous experience using E-learning, once the server gets busy, the entire system is slow. Using Wikiversity and MediaWiki, the traffic of the server never seems to overwhelm the system and cause slow loading and sluggish performance. This increases the accessibility and reliability of MediaWiki. E-learning's information is only available for a semester and then whatever is posted is taken off. On Wikiversity, the information available on it is essentially ever-lasting so long as the user and Wikiversity continues to exist. Using E-learning, the user interface of the system was more cluttered and unorganized than Wikiversity. I believe that Wikiversity has a more user-friendly system. Learning to use MediaWiki was not difficult although typing in the mathematical formulas were very tedious. As for listing MediaWiki experience on a resume, this is a very minor thing. The MediaWiki software is not difficult to use. MATLAB experience, C++ programming experience or web page design experience are far better things to have and list on a resume. Also, a positive aspect of the wiki system is that anyone can edit the information provided by the wiki pages. This allows for universal editing, and accessing. MediaWiki allows for text, images, equations, video, and audio to be displayed on one page all at once, available for immediate use. E-Learning presents the media in file or link form. Text, equations and some images are commonly presented as Word documents on E-Learning. Video, audio and some images are presented as separate files, or as links to those files. Both E-Learning and MediaWiki can present any media type; they just differ on how it is presented. E-Learning has its advantages by presenting the material as files. Printing documents is simpler when they are already in document form, as opposed to trying to print a webpage. Also, if any of the video or audio files were wanted to be on an mp3 player or some other portable device, they can be moved quickly and simply with E-Learning’s file format.

Recommended Software

One helpful Firefox plugin used to edit Wiki articles is wikEd. It features a rich-text WYSIWYG editor with many powerful tools.

Other helpful links included an HTML-to-Wiki converter and an Excel-to-Wiki converter.

For many of the images the open source software, Inkscape, and the vector-based Adobe Illustrator (can be accessed in most UF Computer Labs) were used. And for the equations the latex equation editor proved very helpful.

Contributing Members

Jin Yu Guan Eas4200c.f08.aeris.guan 23:28, 8 December 2008 (UTC)

Delvis Stoute Eas4200c.f08.spars.stoute 06:25, 9 December 2008 (UTC)

Nelson Caceres Eas4200c.f08.aeris.caceres 07:59, 9 December 2008 (UTC)

Josh Holladay Eas4200c.f08.aeris.holladay 16:24, 9 December 2008 (UTC)

Anett Krammer Eas4200c.f08.aeris.krammer 18:27, 9 December 2008 (UTC)

Jesse McRae Eas4200c.f08.aeris.mcrae 02:50, 10 December 2008 (UTC)