User:Eas4200c.f08.blue.a/Homework3

Aerospace Structures Team Blue Homework Report # 3

Lecture Notes
Assigned Reading from Mechanics of Aircraft Structures (C. T. Sun): pg. 451 - Shear Panel pg 85-86 - Shear Flow (Eq. 3.49, 3.50) pg. 20 - Fig. 2.2 pg. 64 - Fig, 3.2 pg 115 - Fig. 4.1

The lecture began with a discussion of shear in non-uniform objects such as the one in Figure 1 below.



Figure 1

Then, curved panels were discussed and reference was made to equations 1.4 and 1.5 (Sun, 5). To illustrate the concept of shear in a curved panel, problem 1.5 from the text book (Sun, 17) was discussed.

Problem 1.5
Problem statement: Derive the relations given by equations 1.4 $$\left( V_{x}=\tau ta\right)$$ and 1.5 $$\left( V_{y}=\tau tb\right)$$ in a curved panel as in Figure 2 below.



Figure 2

A Closer look is taken at the z,y plane and the shear flow on the edge of the curve as seen below:



Figure 3

Now let q be equal to the shear flow traveling from point A to point B as shown by the direction of the arrows above in Figure 3.

$$q= \left( \tau \right) \left( t \right)$$

The infinitesimal change of shear force $$\left( d\vec{F}\right)$$ is shown below as

$$d\vec{F}=qd\vec{l}=q(dl_{y}\hat{j}+dl_{z}\hat{k})=q(dlcos\theta \hat{j}+dlsin\theta \hat{k})=qdydz$$

A closer look at $$d\vec{F}$$is shown in Figure 4 below.



Figure 4

The Resultant shear force vector is found as

$$\vec{F}=\int_{B}^{A}{d\vec{F}}=q\left(\left(\int_{B}^{A}{dy}\right)\hat{j}+\left(\int_{A}^{B}{dz} \right) \hat{k}   \right)$$

Substituting equations 1.4 and 1.5,

$$\vec{F}=q\left(a\hat{j} +b\hat{k}\right)=F_{y}\hat{j}+F_{z}\hat{k}$$

where

$$F_{y}=\left( q\right)\left(a\right)$$ as in equation 3.49b

and

$$F_{z}=\left( q\right)\left(b\right)$$ as in equation 3.49c

yielding

$$\frac{F_{y}}{F_{z}}=\frac{a}{b}$$

The Resultant magnitude is

$$\left|\left|\vec{F} \right| \right|=\sqrt{\left[\left(F_{y} \right)^{2} +\left(F_{z} \right)^{2}\right]}=q\left[a^ {2} +b^{2}\right]^{1/2}$$

$$R=\left|\left|\vec{F} \right| \right|$$

$$T=2\left( q\right)\left( A\right)$$

Now, a thin walled cross section can be analyzed in the same way as the curved panel. The thin walled cross section is seen in the z,y plane below.



Figure 5

a close up of the section in question:



Figure 6

It is assumed that the Torque (T) is in the x direction ie. $$\vec{T}=T\hat{i}$$

$$\vec{T}=\vec{r}\bowtie d\vec{F}$$

It follows that:

$$dT = \rho dF=\rho \left(qdl \right)$$

and therefore:

$$T=\oint_{}^{}{dT}=q\oint_{}^{}{\rho dl}$$

From Figure 6 above, it can be seen that $$   dA=\frac{1}{2}\rho dl$$

Therefore:

$$T=\oint_{}^{}{dT}=q\oint_{}^{}{2dA}=2q\oint_{}^{}{dA}=2q\bar{A}$$

Where $$\bar{A}$$ is equal to the average area as measured from the middle of the thickness dl.


 * {| class="toccolours collapsible collapsed" width="100%" style="text-align:left"

! Homework Solution: Triangle Area Proof

Problem: prove that $$A_{BEC}=\frac{1}{2}(ED)(BC)$$

Solution: The area of the triangle BED is calculated from the formula for the area of a right triangle: $$A=\frac{1}{2}bh$$ b= lenght of horizontal base h= length of vertical height Thus, the area of triangle BED is:

$$A_{BED}=\frac{1}{2}(ED)(BD)$$

and the area of the interior right triangle CED is:

$$A_{CED}=\frac{1}{2}(ED)(CD)$$

Therefore, the area of triangle BEC can be found by subtracting these two areas:

$$A_{BEC}=\frac{1}{2}(ED)(BD)-frac{1}{2}(ED)(CD)=\frac{1}{2}(ED)(BC)$$ |}

For the open thin walled cross section:



This means that: $$T=2q\bar{A}$$ and that the reation: J=\frac{1}{2}\pi a^{4}

Uniform Bar with Circular Cross Section
Nonwarping case: Cross Section behaves as rigid disks



For this case, Torque is calculated thus:

$$T=\int_{A}^{}{\int_{}^{}{}r\tau dA}$$  with   $$\tau =G\gamma $$

and define:

$$\gamma =\frac{rd\alpha }{dx}$$

Therefore: $$T=\int_{A}^{}{rG(r\theta )dA}$$ where G and $$\theta$$ are independent of variables y and z.

This means that the above formula can be written as: $$T=G\theta (\int_{A}^{}{r^{2}dA})$$.

The integral enclosed in brackets is known as the polar moment of inertia, J. if the radius of a circular cross section is defined as a, then J can written as:

$$J=\frac{1}{2}\pi a^{4}$$

Now consider the case of a hollow circular cross section with thin walls (t<<a)

Let: $$r_{inner}=a$$ and $$r_{outer}=b$$

This yields a polar moment of inertia of:

$$J=\frac{1}{2}\pi (b^{4}-a^{4})$$

and with some manipulation:

$$J=\frac{1}{2}\pi(b-a)(b+a) (b^{2}-a^{2})$$

Where (b-a) can be approximated as t and (b+a) is approximately $$2r\bar{}$$.

in addition we assume:

$$b^{2}\cong r\bar{}^{2}$$ $$a^{2}\cong r\bar{}^{2}$$

Substituting yields the eqution:

$$J=2\pi tr\bar{}^{3}$$

rearranging gives the equation:

$$J=2\pi ^{-\frac{1}{2}}t(\pi r\bar{}^{2})^{\frac{3}{2}}$$

where $$\pi r\bar{}^{2}=A\bar{}$$. Therefore, J is proportional to $$A\bar{}^{\frac{3}{2}}$$ with $$(2\pi ^{-\frac{1}{2}}t)$$ as the proportionality factor.

Open Thin Walled Cross Sections
Imagine a uniform bar with circular cross section that is said to be non-warping:





The cross section behaves as a rigid disk

The relation we will take as the connection between Torque that is applied and the resulting twist produced with respect to the strength of the material is called Hooke's Law. Hooke's Law can be derived from the equation of torque by setting certain variables as constants.



where



and



which can then be represented as



or



where a=the radius of the circular cross section

You will notice this equation bears a striking resemblance to the equation for the 2nd area moment of inertia of a body. In this case, though, it is in polar coordinates, and is considered the 2nd polar area moment of inertia. The 2nd polar area moment of inertia is also denoted as J.

Hollow Thin Wall Cross Section
Ri = inner radius | Ro = outer radius





are approximations which give us



or

by arranging the equation this way, we can now see that J is proportional to the area to the 3/2 power.

$$\theta = \frac{\alpha}{x} $$ rate of twist    (1)

$$u_{y}=-\theta \times z $$ is the horizontal portion and

$$u_{z}=+(PP' )\cos \beta$$ is the vertical portion

$$=\alpha y_{p}$$ where $$\alpha = \theta x$$

$$u_{z}=+\theta \simeq y$$    (2)

z is the displacement of a point on the cross section

For the warping displacement along the x-axis is

$$ u_{x} = \theta \psi(y, z)$$    (3)

The kinematic assumptions are equations (1), (2), and (3)

After the following road map, this derivation will continue.

Road map for torsional analysis of an aircraft wing
A. Kinematic Assumption

B. Strain Displacement Relationship

C. Equilibrium Equations

D. Prandth Stress Function $$/phi$$

E. Strain Compatibility Equations

F. Equation for $$\phi$$

G. Boundary Condition for $$\phi$$

H. $$T = 2\int \int \phi dA $$

T = GJ$$\theta$$ J = moment of inertia for a circular cross section Moment of inertia for a non circular cross section is    $$J = \frac{-4}{\bigtriangledown ^{2}\phi}\int \int \phi dA$$

I. Thin walled cross section

Ad-Hoc Assumption Formal Derivation $$T = 2q\bar{A}$$

J. Twist angle $$\theta$$: Method 1

$$\theta = \frac{1}{2G\bar{A}} \oint_{}^{}{}\frac{q}{t}ds$$ s = curvilinear coordinate along this wall

K. Sec. 3.6 on multi cell thin walled cross-section

Continuing the discussion on multicell torsion: The two important equations for single cell torsion are $$T=2*q*A$$ and $$\theta=\frac{1}{2(G)(A)}\oint\frac{q}{t}dl$$

Which is the torque and angle of twist respectively of the cell.

When moving to a multicell view in which:

cell $$i=1,....,n_{cell}$$

Torque on the multicell system now becomes equal to the sum of all the individual cell torques.

$$T=2*\sum{q_i*A_i}$$

in which:

$$q_i$$ is the shear flow in cell i

$$A_i$$ is the "average" area in cell i

so therefore:

$$T_i=2*q_i*A_i$$

so the torque on the for the entire multicell system would then be the sum of all the $$T_i$$

$$T=\sum_{n=1}^{n_c}T_i$$

the same approach is used for the angle of twist

$$\theta$$

in which:

$$\theta= \theta_1,.....,\theta_i$$

Therefore:

$$\theta_i=\frac{1}{2*G_i*A_i}\oint\frac{q_i}{t_i}dl$$ in which $$G_i$$ is the shear modulus of cell i   and $$t_i$$ is the thickness of cell i

So therefore the angl$$e$$ of twist for the entire multicell system would then be the sum of all the

$$\theta_i$$ $$\theta=\sum_{n=1}^{n_c}\theta_i$$

Example

$$t_1=.008m$$ $$t_2=t_3=.01m$$ $$a=4m$$ $$b=2m$$

Therefore the area will equal the sum of the semicircle and triangle which make up the overall airfoil shape.

$$A=\frac{\pi*(\frac{b}{a})^2}{2}+\frac{1}{2}*b*a=5.5708m $$ and $$T=2*q*A$$

So therefore:

$$q=\frac{T}{2*A}$$

in which T is the variable, so the twist angle will equal:

$$\theta=\frac{1}{2*G*A}\sum_{j=1}^{3}\frac{q_j*l_j}{t_j}$$

in which l is the length of the path.

NACA Airfoil Matlab Code
functi on NACA_Airfoil_Generator(m,p,t)

% This function takes the 4 digit NACA airfoil digits and plots the % airfoil. It then plots the centroid of the airfoil and calculates the % Area of the airfoil.

% m = maximum camber in percentage of chord length. % p = position of maximum camber in tenths of the chord length. % t = maximum airfoil thickness in percentage of the chord.

% Chord Length = 1

% Converting the NACA numbers into useable values for calculations.

m = 0.01*m; p = 0.1*p; t = 0.01*t;

% Creating the x input arrays.

x = [0:0.01:1];

x_c1 = [0:0.01:p]; x_c2 = [p+0.01:0.01:1];

% Calculating the y coordinates of the mean camber line.

y_c1 = (m/(p^2)).*(2*p.*x_c1-(x_c1.^2)); y_c2 = (m/((1-p)^2)).*((1-2*p)+(2*p.*x_c2)-(x_c2.^2));

y_c = cat(2, y_c1, y_c2)

% Computing the thickness distribution of the airfoil

y_t = 5*t*((0.2969.*sqrt(x))-0.1260.*x-(0.3516.*(x.^2))+(0.2843.*(x.^3))-(0.1015.*(x.^4)))

% Taking the derivative of the camber line. dy_c1 = (m/(p^2)).*(2*p-(2.*x_c1)); dy_c2 = (m/((1-p)^2)).*((1-2*p)+(2*p)-(2.*x_c2));

dy_c = cat(2,dy_c1, dy_c2)

% Computing the angle between the x-axis to the tangent of the mean camber % line.

theta = atan(dy_c)

% Computing the upper and lower coordinates of the airfoil

x_u = x - y_t.*sin(theta); y_u = y_c + y_t.*cos(theta);

x_l = x + y_t.*sin(theta); y_l = y_c - y_t.*cos(theta);

% Plotting the mean camber line, upper coordinates, and lower coordinates % of the airfoil. Calculating and plotting the centroid as a cross-hair.

% Calculating the centroid ns = 100; x_centroid = 0; y_centroid = 0; total = 0;

% Upper airfoil section for i = 1:ns dt = (x_u(i+1)-x_u(i))*(y_u(i+1)+y_u(i))*0.5; d_centroid_x = (x_u(i+1)+x_u(i))*0.5; d_centroid_y = (y_u(i+1)+y_u(i))*0.5; x_centroid = x_centroid + (dt*d_centroid_x); y_centroid = y_centroid + (dt*d_centroid_y); total = total + dt; end

% Lower Airfoil Section for i = 1:ns dt = (x_l(i+1)-x_l(i))*(y_l(i+1)+y_l(i))*0.5; d_centroid_x = (x_l(i+1)+x_l(i))*0.5; d_centroid_y = (y_l(i+1)+y_l(i))*0.5; x_centroid = x_centroid + (dt*d_centroid_x); y_centroid = y_centroid + (dt*d_centroid_y); total = total + dt; end

centroid_x = x_centroid/total centroid_y = y_centroid/total

% Plotting

plot(x,y_c) hold on plot(x_u,y_u) plot(x_l,y_l) plot(centroid_x, centroid_y, '+') axis on hold off

% Calculating the value of A-bar.

ns = 100; x_cent = 0; y_cent = 0; a_bar = 0;

% Upper airfoil section for i = 1:ns r = [x_u(i) - x_cent, y_u(i) - y_cent, 0]; dl = [x_u(i+1) - x_u(i), y_u(i+1) - y_u(i), 0]; d_abar = 0.5.*cross(dl,r); a_bar = a_bar+d_abar; end

% Now for the lower airfoil section for i = 1:ns r = [x_l(i) - x_cent, y_l(i) - y_cent, 0]; dl = [x_l(i+1) - x_l(i), y_l(i+1) - y_l(i), 0]; d_abar = 0.5.*cross(dl,r); a_bar = a_bar+d_abar; end

From the code, the calculated A-bar for the NACA 2415 airfoil is 0.1020.

Team Members that Contributed to this Report
Mark Barry Eas4200c.f08.blue.a 04:09, 8 October 2008 (UTC)

Jeffrey Shea Eas4200c.f08.Blue.E 17:38, 7 October 2008

Christian Garabaya Eas4200c.f08.blue.d 12:58, 8 October 2008 (UTC)

Thomas Sexton EAS4200c.f08.blue.c 13:14, 8 October 2008 (UTC)

Marc Pelletier EAS4200c.f08.blue.b 16:16, 8 October 2008 (UTC)

Clay Robertson eas4200c.f08.blue.f 23:25, 8 October 2008 (UTC)