User:EAS4200C.F08.WIKI.A/Homework4

=Quadrature Method=

Quadratures can be used to find the average area of foils. Below is an image of the quadrature of a foil.



$$\overline{A} = \overline{A}_1 + \overline{A}_2 = \overline{A}_1 - \begin{vmatrix}\overline{A}_2\end{vmatrix}$$
 * $$ \overline{A}_1 > 0 $$
 * $$ \overline{A}_2 > 0 $$




 * $$\overrightarrow{r}x\overrightarrow{PQ}\; goes\; in\; negative\; x\; direction\; $$

Using Trapezoidals to Integrate
There are problems with using this method because of the curvature of the foil. When you split up the area of the foil, not every section will end up looking like a trapezoid. This method should not be used for foils.



Shear Flow in a Single Cell Airfoil
Below is an image of a single cell airfoil.



The shear flow is constant, therefore:


 * $$ q=q_1=q_2=q_3$$

Rate of Twist Angle θ

$$\theta = \frac{1}{2G\overline{A}}\; q\sum_{j=1}^{}{\frac{l_j}{t_j}}$$

$$\theta = \frac{1}{2G\overline{A}}\; q[\frac{\pi (\frac{b}{2})}{t_1} + \frac{a}{t_2} + \frac{\sqrt{a^2 + b^2}}{t_2}]$$

$$\theta = \frac{1}{2G\overline{A}}\; q[\frac{\pi bt_2t_3 + 2at_1t_3 + 2t_1t_2\sqrt{a^2 + b^2}}{t_1t_2t_3}]$$

Shear Stress
Max Shear Stress τmax

$$\tau_{max} = \frac{q}{min(t_1, t_2, t_3)}$$


 * $$if\; \tau_{max} = \tau_{all}$$


 * $$and\; since\; q = \frac{T}{2\overline{A}} \; then$$

$$T_{all} = 2\bar{A}\tau_{all}\;min(t_1, t_2, t_3)$$

=Multicell Airfoil=

In this section we will examin a multicel airfoil. We are interested to find the rate of twist as a function of the torque $$\displaystyle T$$ and the torsional constant $$\displaystyle J$$. The figure depicted below is the object that will be in concideration.



We first have to express the torque of the while object. To do so we will can use the principal of superposition. In other words the torque of the whole object is the summation of the torques from each cell, $$\displaystyle C_1$$ and $$\displaystyle C_2$$.

$$\displaystyle \sum_{i=1}^{N_c}{2q_i\bar{A_i}}= 2q_1\bar{A_1} + 2q_2\bar{A_2}$$

We now must find the rate of twist for each cell. The formulation for cell one is shown below.

$$\displaystyle \theta_1 = \frac{1}{2G\bar{A_1}}\oint_{a}^{a}{\frac{q_1}{t_1} ds} = \frac{1}{2G\bar{A_1}}[\frac{2q_1a}{t_1} + \frac{q_1c}{t_1} + \frac{(q_1-q_2)c}{t_{12}}]$$

When all known variables are put in and noting that $$\displaystyle \bar{A_1} = ac = (30cm)(40cm) = 1200cm^2$$ equals. the expression for the rate of twist for cell one then becomes;

$$\theta_1 = \frac{1}{Gcm^2}[.181q_1 - .042q_2]$$

The same is done for cell two and is performed below.

$$\displaystyle \theta_2 = \frac{1}{2G\bar{A_2}}\oint_{a}^{a}{\frac{q_2}{t_2} ds} = \frac{1}{2G\bar{A_2}}[\frac{2q_2b}{t_2} + \frac{q_2c}{t_2} + \frac{(q_2-q_1)c}{t_{12}}]$$

Substituting values and noting that $$\displaystyle \bar{A_2} = bc = (60cm)(40cm) = 2400cm^2$$, the expression reduces to;

$$\theta_2 = \frac{1}{Gcm^2}[-.021q_1 + .062q_2]$$

Now it is imperative to notice that the rates of twist for each cell are equivalent. Setting them equal to each other and solving for either $$\displaystyle q_1$$ or $$\displaystyle q_2$$ is the next step in solving this problem. the result is as follows. We chose to solve for $$\displaystyle q_1$$.

$$ \frac{1}{Gcm^2}[.181q_1 - .042q_2] = \frac{1}{Gcm^2}[-.021q_1 + .062q_2]$$  $$\displaystyle q_1= .65 q_2$$

Now we pace this expression for $$\displaystyle q_1$$ into the torque equation that was first derived in order to obtain the rate of twist for cell two $$\displaystyle q_2$$ in terms of the torque $$\displaystyle T$$. It is then possible to obtain the rate of twist for cell one $$\displaystyle q_1$$ in terms of the torque.

$$\displaystyle T = 2q_2(.65)(1200cm^2) + 2q_2(2400cm^2) $$  $$\displaystyle q_2 = \frac{1.57x10^{-4}T}{cm^2}$$  $$\displaystyle q_1 = \frac{1.02x10^{-4}T}{cm^2}$$

Now we must put these expressions back into either one of the rates of twist equations. Doing so will put the rate of twist in terms of the torque. This is done below.

$$\theta_1 = \frac{1}{Gcm^2}[.181q_1 - .042q_2] = \frac{1}{Gcm^2}[.181(\frac{1.02x10^{-4}T}{cm^2}) - .042(\frac{1.57x10^{-4}T}{cm^2})]$$  <p style="text-align:center;">$$\theta_1 = \frac{1.18x10^{-5}T}{Gcm^4} = \theta_2$$

Now we must relate the angle of twist to the torsional constant. We know that the angle of twist is equal to the following.

<p style="text-align:center;">$$\displaystyle \theta = \frac{T}{GJ}$$

We can conclude that

<p style="text-align:center;"> $$\displaystyle \frac{1}{J} = \frac{1.18x10^{-6}}{cm^4}$$

And that

<p style="text-align:center;">$$\displaystyle J = \frac{1}{1.18x10^{-6}} = 84745.8cm^4$$

=Ad-Hoc Derivation=

The engineering, or ad-hoc, derivation of the rate of twist of the following equation

<p style="text-align:center;"> $$\displaystyle \theta = {1 \over 2G \bar A} \oint {q \over t} ds $$

applies to a uniform bar with a non-circular cross section.



The displacement of segment $$\displaystyle PP' $$ due to an angle $$\displaystyle \alpha $$ is defined as

<p style="text-align:center;"> $$\displaystyle {PP' \over OP} = \tan \alpha $$

But for a small angle $$\displaystyle \alpha $$, the following approximation is appropriate:

<p style="text-align:center;"> $$\displaystyle \tan \alpha \approx \alpha $$

Hence,

<p style="text-align:center;"> $$\displaystyle {PP' \over OP} \approx \alpha $$

The projected displacement of $$\displaystyle PP' $$ in the direction perpendicular to $$\displaystyle OP' $$ can be written

<p style="text-align:center;"> $$\displaystyle PP'' = PP' \cos \alpha $$

By substitution

<p style="text-align:center;"> $$\displaystyle PP'' = (OP \tan \alpha) \cos \alpha $$

By the associative property of multiplication

<p style="text-align:center;"> $$\displaystyle PP'' = (OP \cos \alpha) \tan \alpha $$

Recall that $$\displaystyle OP = r $$ and $$\displaystyle OP'' = \rho $$, thus we can express it as

<p style="text-align:center;"> $$\displaystyle PP'' = (r \cos \alpha) \tan \alpha $$

Strain, $$\displaystyle \gamma $$ can be written as

<p style="text-align:center;"> $$\displaystyle \gamma = {PP' \over dx} = {\rho \alpha \over dx} = \rho \theta $$

where, as previously stated, $$\displaystyle \theta $$ is the rate of twist.

The expression:

<p style="text-align:center;">$$\theta = \frac{\partial \alpha }{\partial x}$$

can be used to express the rate of twist of a cross section.

Hooke's law can be represented as:

<p style="text-align:center;">$$\tau = G\gamma = G\rho \theta$$

<p style="text-align:center;"> $$\tau (s)= G\gamma = G\rho(s) \theta (x)$$

With variables as defined in the previous diagram.

Integrate along $$\textit{C}$$

<p style="text-align:center;">$$ \oint_\textit{C}{}^{}{\tau (s)ds} = G\theta (x)\oint_{\textit{C}}^{}{\rho (s)ds}$$

where

<p style="text-align:center;">$$ \tau (s)ds = \frac{q(s)}{t(s)}$$

$$\rho (s) = 2\bar{A}$$

To get an expression for $$\theta$$ as discussed in the previous lecture.

The question, 'What is ad-hoc about the above derivation of $$\theta$$ expression on p 20.2 and the derivation of $$T = 2q\bar{A}$$?' was posed in the discussion of the equation. The answer can be provided in two parts:

1) Strain ($$\gamma$$) must be obtained using the displacement of $$P$$ in the direction tangent to $$\textit{C}$$ at $$P$$, but $$PP''$$ in the previous diagram is not necessarily tangent to $$\textit{C}$$ (but is actually close, due to the small angles involved).

2) $$\tau = \frac{q}{t}$$ obtained from ad-hoc assumption that $$\tau$$ was uniform across wall thickness, this is assuming that the wall thickness is very thin, or else this assumption breaks down.

Now formal justification (derivation) of elasticity theory.

This section refers back to the road map that was previously discussed.

A. Kinematic Assumption

<p style="text-align:center;">$$u_{x}(y, z) = \theta \Psi (y,z)$$

Where $$\theta$$ is considered constant with respect to x (due to it being a uniform bar).

<p style="text-align:center;">$$u_{y}(x, z) = -\theta xz$$ <p style="text-align:center;"> $$u_{z}(x, y) = \theta xy$$

To transform equations in Sun[2006] to those using our unified notation, we use cyclic permutation

<p style="text-align:center;"> $$\varepsilon _{xx} = \varepsilon _{yy} = \varepsilon _{zz} = \gamma _{yz} = 0$$

$$\varepsilon _{xx} = \frac{\partial u_{x}}{\partial x}(y,z) = 0$$

$$\gamma _{yz} = \frac{\partial u_{y}}{\partial z}(x,z) + \frac{\partial u_{z}}{\partial y}(x,y) = -\theta x + \theta x = 0$$ HW:

<p style="text-align:center;">$$\varepsilon _{yy} = \frac{\partial u_{y}}{\partial y}(x,z) = 0$$

$$\varepsilon _{zz} = \frac{\partial u_{z}}{\partial z}(x,y) = 0$$

These results are obtained due to the fact that both derivatives are taken with respect to a variable that does not show up in either function being differentiated.

=Strain=

Normal strains

$$\varepsilon _{xx}=\frac{\sigma _{xx}}{E}-\frac{\nu \sigma _{yy}}{E}-\frac{\nu\sigma _{zz} }{E} $$

$$\varepsilon _{yy}=-\frac{\nu\sigma _{xx}}{E}+\frac{ \sigma _{yy}}{E}-\frac{\nu\sigma _{zz} }{E} $$

$$\varepsilon _{zz}=-\frac{\nu\sigma _{xx}}{E}-\frac{\nu \sigma _{yy}}{E}+\frac{\sigma _{zz} }{E} $$

Shear Strains

$$\gamma _{xy} =2 \zeta _{xy}= \frac{\tau_{xy}}{G}=\frac{\sigma_{xy}}{G}$$

$$\gamma _{yz} =2 \zeta _{yz}= \frac{\tau_{yz}}{G}=\frac{\sigma_{yz}}{G}$$

$$\gamma _{zx} =2 \zeta _{zx}= \frac{\tau_{zx}}{G}=\frac{\sigma_{zx}}{G}$$

Unroll this matrix along the diagonal and then up and around counterclockwise.

$$\epsilon = \begin{vmatrix} \varepsilon_ {11} & \varepsilon_{12} & \varepsilon _{13}\\ \varepsilon _{21} &\varepsilon _{22} &\varepsilon _{23} \\ \varepsilon _{31} &\varepsilon _{32} &\varepsilon _{33} \end{vmatrix}$$

Due to the symmetry, you may rewrite this unrolled matrix into a $$6_x 1$$ with a Voigt Notation:

$$ \left\{\varepsilon _{ij} \right\}_{6x1}= \begin{Bmatrix} \varepsilon _{11}\\ \varepsilon _{22}\\ \varepsilon _{33}\\ \varepsilon _{23}\\ \varepsilon _{13}\\ \varepsilon _{12} \end{Bmatrix}$$

$$\left\{\sigma _{ij} \right\}_{6x1}= \begin{Bmatrix} \sigma _{11}\\ \sigma _{22}\\ \sigma _{33}\\ \sigma _{23}\\ \sigma _{13}\\ \sigma _{12} \end{Bmatrix}$$

Hooke's Law for isotropic elasticity:

$$\begin{Bmatrix} \varepsilon _{11}\\ \varepsilon _{22}\\ \varepsilon _{33}\\ \gamma _{23}\\ \gamma _{13}\\ \gamma _{12} \end{Bmatrix}=

\begin{vmatrix} \frac{1}{E}& \frac{-\nu }{E} & \frac{-\nu }{E} & 0 &0 &0 \\ \frac{-\nu }{E}& \frac{1}{E} &\frac{-\nu }{E}  & 0 & 0 & 0\\ \frac{-\nu }{E} & \frac{-\nu }{E} & \frac{1}{E} & 0 &0  &0 \\ 0& 0 & 0 & \frac{1}{G} & 0 & 0\\ 0& 0 & 0 & 0 & \frac{1}{G} &0 \\ 0& 0 & 0 & 0 & 0 & \frac{1}{G} \end{vmatrix} \begin{Bmatrix} \sigma _{11}\\ \sigma _{22}\\ \sigma _{33}\\ \sigma _{23}\\ \sigma _{13}\\ \sigma _{12} \end{Bmatrix}$$

It can also be written as

$$\begin{Bmatrix} \varepsilon _{11}\\ \varepsilon _{22}\\ \varepsilon _{33}\\ \varepsilon _{23}\\ \varepsilon _{13}\\ \varepsilon _{12} \end{Bmatrix}=

\begin{vmatrix} \frac{1}{E}& \frac{-\nu }{E} & \frac{-\nu }{E} & 0 &0 &0 \\ \frac{-\nu }{E}& \frac{1}{E} &\frac{-\nu }{E}  & 0 & 0 & 0\\ \frac{-\nu }{E} & \frac{-\nu }{E} & \frac{1}{E} & 0 &0  &0 \\ 0& 0 & 0 & \frac{1}{2G} & 0 & 0\\ 0& 0 & 0 & 0 & \frac{1}{2G} &0 \\ 0& 0 & 0 & 0 & 0 & \frac{1}{2G} \end{vmatrix} \begin{Bmatrix} \sigma _{11}\\ \sigma _{22}\\ \sigma _{33}\\ \sigma _{23}\\ \sigma _{13}\\ \sigma _{12} \end{Bmatrix}$$

Poisson's ratio values for different materials
Poisson's Ratio table []

=MATLAB=

Idealized Single Cell Airfoil
function IdealizedSingleCell

a = 4; b = 2; r = b/2; t1 = .008; t2 = .01; t3 = .01;

y = linspace(0, -r, 25);

y(26:50) = linspace(-r, 0, 25);

y(51) = a; y(52) = 0;

z = sqrt(r^2 - y(1:25).^2);

z(26:50) = -sqrt(r^2 - y(26:50).^2);

z(51) = -r; z(52) = r;

plot(y,z)

A_bar = pi/2*(b/2)^2 + .5*b*a;

ThetaG = 1/(2*A_bar)*(pi*(b/2)/t1 + a/t2 + sqrt(a^2+b^2)/t3); % *q

J = 2*A_bar/ThetaG

EDU>> IdealizedSingleCell

J =

0.1001



Rectangular 2-Cell Airfoil
function Rectangle2Cell

a = .3; b = .6; c = .4;

t1 = .003; t2 = .005; t12 = .004;

A1_bar = a*c; A2_bar = b*c;

X = 2*a/t1 + c/t1 + c/t12;

Y = 2*b/t2 + c/t2 + c/t12;

Z = -c/t12;

alpha = (A2_bar*X - A1_bar*Z)/(A1_bar*Y - A2_bar*Z);

J = (4*A1_bar^2 + 4*alpha*A1_bar*A2_bar)/(X + Z*alpha)

EDU>> Rectangle2Cell

J =

8.5507e-004

Circle Validation
function Circle(c,P_0)

R = c/2; check = 1; j = 1;

while check == 1 ns = 2*j; y = linspace(0,c,ns); for i = 1:ns z_u(i) = sqrt(R^2-(y(i)-R)^2); z_l(i) = -sqrt(R^2-(y(i)-R)^2); end

plot(y,z_u,y,z_l)

A_bar = 0; y_bar = 0; z_bar = 0; s = 0; y_u = y;   y_l = y;

for i = 1:ns-1 PQ = [0 y(ns-i)-y(ns-i+1) z_u(ns-i)-z_u(ns-i+1)]; r = [0 y(ns-i+1)-P_0(2) z_u(ns-i+1)-P_0(3)]; a = cross(r,PQ)/2; A_bar = A_bar + a;       y_bar = y_bar + (y_u(i)+y_u(i+1))/2*norm([y_u(i+1) z_u(i+1)]-[y_u(i) z_u(i)]); z_bar = z_bar + (z_u(i)+z_u(i+1))/2*norm([y_u(i+1) z_u(i+1)]-[y_u(i) z_u(i)]); s = s + norm([y_u(i+1) z_u(i+1)]-[y_u(i) z_u(i)]); end

for i = 1:ns-1 PQ = [0 y(i+1)-y(i) z_l(i+1)-z_l(i)]; r = [0 y(i)-P_0(2) z_l(i)-P_0(3)]; a = cross(r,PQ)/2; A_bar = A_bar + a;       y_bar = y_bar + (y_l(i)+y_l(i+1))/2*norm([y_l(i+1) z_l(i+1)]-[y_l(i) z_l(i)]); z_bar = z_bar + (z_l(i)+z_l(i+1))/2*norm([y_l(i+1) z_l(i+1)]-[y_l(i) z_l(i)]); s = s + norm([y_u(i+1) z_u(i+1)]-[y_u(i) z_u(i)]);

end

NS(j) = ns; A_BAR(j) = A_bar(1);

if abs(pi*R^2-A_bar(1))/(pi*R^2)<.01 check = 0; else j = j+1; end end

y_bar = y_bar/s; z_bar = z_bar/s;

A_bar = A_BAR(j)

subplot(2,1,1), plot(y_u,z_u,y_l,z_l,y_bar,z_bar,'x'), xlabel('y'), ylabel('z'), title('NACA 2415 Airfoil')

subplot(2,1,2), plot(NS,A_BAR), xlabel('ns'), ylabel('A_b_a_r'), title('ns vs. A_b_a_r')

EDU>> Circle(1,[0 0 0])

A_bar =

0.7779



NACA 2415 Single Cell Airfoil
function Airfoil(NACA,c,P_0)

m = floor(NACA/1000)/100;

p = (floor(NACA/100)-floor(NACA/1000)*10)/10;

ta = (NACA-floor(NACA/100)*100)/100;

ts = .002;

check = 1; j = 1;

while check ==1 ns = 5*j; y = linspace(0,1,ns);

for i = 1:ns if y(i)<=p z_c(i) = (m/p^2)*(2*p*y(i)-y(i)^2); theta(i) = atan((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); theta(i) = atan((m/(1-p)^2)*(2*p-2*y(i))); end z_t(i) = 5*ta*(.2969*y(i)^.5-.1260*y(i)-.3516*y(i)^2+.2843*y(i)^3-.1015*y(i)^4);

y_u(i) = (y(i) - z_t(i)*sin(theta(i)))*c; z_u(i) = (z_c(i) + z_t(i)*cos(theta(i)))*c;

y_l(i) = (y(i) + z_t(i)*sin(theta(i)))*c; z_l(i) = (z_c(i) - z_t(i)*cos(theta(i)))*c; end

Area = 0; y_bar = 0; z_bar = 0; s = 0;

for i = 1:ns-1 PQ = [0 y_u(ns-i)-y_u(ns-i+1) z_u(ns-i)-z_u(ns-i+1)]; r = [0 y_u(ns-i+1)-P_0(2) z_u(ns-i+1)-P_0(3)]; a = cross(r,PQ)/2; Area = Area + a;       y_bar = y_bar + (y_u(i)+y_u(i+1))/2*norm([y_u(i+1) z_u(i+1)]-[y_u(i) z_u(i)]); z_bar = z_bar + (z_u(i)+z_u(i+1))/2*norm([y_u(i+1) z_u(i+1)]-[y_u(i) z_u(i)]); s = s + norm([y_u(i+1) z_u(i+1)]-[y_u(i) z_u(i)]); end

for i = 1:ns-1 PQ = [0 y_l(i+1)-y_l(i) z_l(i+1)-z_l(i)]; r = [0 y_l(i)-P_0(2) z_l(i)-P_0(3)]; a = cross(r,PQ)/2; Area = Area + a;       y_bar = y_bar + (y_l(i)+y_l(i+1))/2*norm([y_l(i+1) z_l(i+1)]-[y_l(i) z_l(i)]); z_bar = z_bar + (z_l(i)+z_l(i+1))/2*norm([y_l(i+1) z_l(i+1)]-[y_l(i) z_l(i)]); s = s + norm([y_u(i+1) z_u(i+1)]-[y_u(i) z_u(i)]); end

NS(j) = ns; A_BAR(j) = Area(1); J(j) = 4*Area(1)^2*ts/s;

if j ~= 1 if abs((J(j)-J(j-1))/J(j-1))*100<1 check = 0; end end j = j+1; end

y_bar = y_bar/s; z_bar = z_bar/s;

A_bar = A_BAR(j-1)

subplot(2,1,1), plot(y_u,z_u,y_l,z_l,y_bar,z_bar,'x'), xlabel('y'), ylabel('z'), title('NACA 2415 Airfoil')

subplot(2,1,2), plot(NS,J), xlabel('ns'), ylabel('Torsional Constant, J'), title('ns vs. J')



NACA 2415 3-Cell Airfoil
function Airfoil3Cell(NACA,c,P_0)

m = floor(NACA/1000)/100;

p = (floor(NACA/100)-floor(NACA/1000)*10)/10;

ta = (NACA-floor(NACA/100)*100)/100;

ts = .002; tp = .003;

check = 1; j = 1;

while check ==1 ns = 20*j; y = linspace(0,1,ns);

for i = 1:ns if y(i)<=p z_c(i) = (m/p^2)*(2*p*y(i)-y(i)^2); theta(i) = atan((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); theta(i) = atan((m/(1-p)^2)*(2*p-2*y(i))); end z_t(i) = 5*ta*(.2969*y(i)^.5-.1260*y(i)-.3516*y(i)^2+.2843*y(i)^3-.1015*y(i)^4);

y_u(i) = (y(i) - z_t(i)*sin(theta(i)))*c; z_u(i) = (z_c(i) + z_t(i)*cos(theta(i)))*c;

y_l(i) = (y(i) + z_t(i)*sin(theta(i)))*c; z_l(i) = (z_c(i) - z_t(i)*cos(theta(i)))*c; end

for i = 1:ns if y_u(i) > c/4 P1UpperIndex = i;           break end end

for i = 1:ns if y_l(i) > c/4 P1LowerIndex = i;           break end end

for i = P1UpperIndex:ns if y_u(i) > 3*c/4 P2UpperIndex = i;           break end end

for i = P1LowerIndex:ns if y_l(i) > 3*c/4; P2LowerIndex = i;           break end end

%Cell 1 C1y = (y_u(P1UpperIndex)+y_l(P1LowerIndex))/2; for i = 1:P1UpperIndex C1y(i+1) = y_u(P1UpperIndex-i+1); end for i = 1:P1LowerIndex C1y = [C1y, y_l(i)]; end C1y = [C1y, (y_u(P1UpperIndex)+y_l(P1LowerIndex))/2];

C1z = (z_u(P1UpperIndex)+z_l(P1LowerIndex))/2; for i = 1:P1UpperIndex C1z(i+1) = z_u(P1UpperIndex-i+1); end for i = 1:P1LowerIndex C1z = [C1z, z_l(i)]; end C1z = [C1z, (z_u(P1UpperIndex)+z_l(P1LowerIndex))/2];

%Cell 2 for i = 1:P2UpperIndex-P1UpperIndex C21y(i) = y_u(P2UpperIndex-i+1); end C21y = [C21y, y_l(P1LowerIndex), y_u(P2UpperIndex);

for i = 1:P2UpperIndex-P1UpperIndex C21z(i) = z_u(P2UpperIndex-i+1); end C21z = [C21z, z_l(P1LowerIndex), z_u(P2UpperIndex);

for i=1:P2LowerIndex-P1LowerIndex C22y(i) = y_l(P1LowerIndex+i-1); end C22y = [C22y, y_u(P2UpperIndex), y_l(P1LowerIndex)]

for i=1:P2LowerIndex-P1LowerIndex C22z(i) = z_l(P1LowerIndex+i-1); end C22z = [C22z, z_u(P2UpperIndex), z_l(P1LowerIndex)]

%Cell 3 C3y = (y_u(P2UpperIndex)+y_l(P2LowerIndex))/2; for i = 1:ns-P2LowerIndex C3y(i+1) = y_l(P2LowerIndex+i-1) end for i = 1:ns-P2UpperIndex C3y'(i) = y_u(ns-i+1) end C3y = [C3y, C3y', (y_u(P2UpperIndex)+y_l(P2LowerIndex))/2];

C3z = (z_u(P2UpperIndex)+z_l(P2LowerIndex))/2; for i = 1:ns-P2LowerIndex C3z(i+1) = z_l(P2LowerIndex+i-1) end for i = 1:ns-P2UpperIndex C3z'(i) = z_u(ns-i+1) end C3z = [C3z, C3z', (z_u(P2UpperIndex)+z_l(P2LowerIndex))/2];

%Area & Torsional Constant Area1 = 0; Area21 = 0; Area22 = 0; Area3 = 0;

for i = 1:size(C1y) PQ = [0 C1y(ns-i)-C1y(ns-i+1) C1z(ns-i)-C1z(ns-i+1)]; r = [0 C1y(ns-i+1)-(y_u(P1UpperIndex)+y_l(P1LowerIndex))/2 C1z(ns-i+1)-0]; a = cross(r,PQ)/2; Area1 = Area1 + a;

end

for i = 1:size(C21y) PQ = [0 C21y(i+1)-C21y(i) C21z(i+1)-C21z(i)]; r = [0 C21y(i)-y_l(P1LowerIndex) C21z(i)-z_l(P1LowerIndex)]; a = cross(r,PQ)/2; Area21 = Area21 + a;

end

for i = 1:size(C22y) PQ = [0 C22y(i+1)-C22y(i) C22z(i+1)-C22z(i)]; r = [0 C22y(i)-y_u(P2UpperIndex) C22z(i)-z_u(P2UpperIndex)]; a = cross(r,PQ)/2; Area22 = Area22 + a;

end

for i = 1:size(C3y) PQ = [0 C3y(i+1)-C3y(i) C3z(i+1)-C3z(i)]; r = [0 C3y(i)-(y_u(P2UpperIndex)+y_l(P2LowerIndex))/2 C3z(i)-0]; a = cross(r,PQ)/2; Area3 = Area3 + a;

end

NS(j) = ns; A_BAR(j) = Area(1);

if j ~= 1 if abs((A_BAR(j)-A_BAR(j-1))/A_BAR(j-1))*100<.01 check = 0; end end j = j+1; end

y_bar = y_bar/s; z_bar = z_bar/s;

A_bar = A_BAR(j-1)

P1y = [y_u(P1UpperIndex) y_l(P1LowerIndex)];

P1z = [z_u(P1UpperIndex) z_l(P1LowerIndex)];

P2y = [y_u(P2UpperIndex) y_l(P2LowerIndex)];

P2z = [z_u(P2UpperIndex) z_l(P2LowerIndex)];

subplot(2,1,1), plot(y_u,z_u,y_l,z_l,y_bar,z_bar,'x',P1y,P1z,P2y,P2z), xlabel('y'), ylabel('z'), title('NACA 2415 Airfoil')

subplot(2,1,2), plot(NS,A_BAR), xlabel('ns'), ylabel('A_b_a_r'), title('ns vs. A_b_a_r')



=Individual Contributions=

Christopher Pergola EAS4200C.F08.WIKI.A 20:59, 24 October 2008 (UTC)

Braden Barnes Eas4200c.f08.wiki.b 21:54, 23 October 2008 (UTC)

Chris Fontana Eas4200c.f08.WIKI.E 18:55, 24 October 2008 (UTC)

Eric Viale Eas4200c.f08.wiki.c 19:51, 24 October 2008 (UTC)

Jeff Huenink Eas4200c.f08.wiki.f 20:07, 24 October 2008 (UTC)

Michael Rolle Eas4200c.f08.wiki.d 20:59, 24 October 2008 (UTC)