User:Eas4200c.f08.vqcrew.c/Homework 3

Homework 3

Curved Panels

Figure 1 shows a curved panel with uniform thickness $t$ and constant shear stress $\tau$. This diagram will be analyzed in order to determine the forces acting on the panel in both the y and z directions. The shear flow $q$ is defined as:

Fig 1. Curved Panel

$q=t \tau$

In order to find the total forces acting in the y and z direction, the panel must be divided up into small sections with differential lengths and force vectors. A side view of the panel with dimensions and the differential components is shown in Figure 2. The differential force element, $d\bar{F}$, is equal to the shear flow multiplied by the differential length of the element, or:
$d\bar{F}=qd\bar{L}$

Fig 2. Side view of curved panel

A larger view of the differential length and force element are shown in Figure 3. Figure 3 also shows the components of the differential length element. The differential force vector and length vector can be broken into y and z components and simplified as follows:

Fig 3. Differential force and length units for a section in the curved panel

$d\bar{F}=qd\bar{L}$
$d\bar{F}=q(dL_{y} \hat{j} + dL_{z} \hat{k})$
$d\bar{F}=q(dL cos\theta \hat{j} + dL sin \theta \hat{k})$
$d\bar{F}=q(dy \hat{j} + dz \hat{k})$

By integrating along the path from A to B, we can find the total force acting on the curved panel.

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

The force vector can be broken into y and z components and equated to the above solution to produce equations (1.4) and (1.5):
$\bar{F}=F_{y} \hat{j} + F_{z} \hat{k}$
$\Rightarrow F_{y} = qa = t \tau a$
$\Rightarrow F_{z} = qb = t \tau b$

The relationship between the force components in a curved panel can be expressed as:

$\frac{F_y}{F_z}=\frac{a}{b}$

The magnitude of the force, $\left|\left| \bar{F} \right| \right|$, can be expressed as follows:
$\left|\left| \bar{F} \right| \right|=\sqrt{(F_y)^2+(F_z)^2} = q\sqrt{a^2+ b^2}$

However, $\sqrt{a^2+ b^2}$ is equal to the straight line distance, $d$, from point A to B. Therefore the force vector can also be expressed as:

$\bar{F}=qd$

Torsion in a Closed, Thin-walled Beam

Figure 4 depicts a generic, closed, thin-walled cross-section for a beam under torsion. The coordinate system may be defined from any arbitrary point in the plane of the section.

Fig 4. Closed, thin-walled cross-section of a beam in torsion

The torque vector points in the positive x direction, or:
$\bar{T}=T \hat{i}$

The differential unit of torque can be described as the cross product of the distance from the origin to the differential length segment with the differential force vector, or:
$d \bar{T}= \bar{r} \times d \bar{F}$

The magnitude of this element is:

$dT=\rho dF$

where $\rho$ is the perpendicular distance from the origin to the force vector and $dF=qdL$. To find the total torque, the differential torque element must be integrated about the closed path of the cross-section.

$T=\oint{dT} = \oint{\rho dF}=q\oint{\rho dL}$
However, $\rho dL$ is also equal to twice the area of the triangle swept out by the $r$ vector and $dL$. We can evaluate the integral over the average area, $\bar{A}$, to find the total torque as follows:

$T=2q\int _{\bar{A}}{dA}=2q\bar{A}$

Torsion Acting on a Uniform Bar with Circular Cross-Section

Figure 5 shows a generic case of torsion acting on a solid, circular beam.

This case does not exhibit warping. The only deformation that occurs is contained within the y-z plane. The beam behaves as a series of rigid disks rotating entirely in the plane of the cross-section. Figure 6 displays the reaction of the beam to a given torque.

Fig 5. Torsion acting on a uniform beam with solid circular cross-section.

The general formula for the torque is:

$T = \iint_A r \tau dA$

The shear stress is defined as the product of the shear modulus, G, and the angle of twist, $\gamma$, (Hooke's Law).

$\tau = G \gamma$

$\displaystyle \gamma = \displaystyle \frac{rd\alpha}{dx} = r\theta$

Fig 6. Graphical representation of the effect of torsion on a solid circular beam.

Define $\theta$ as rate of twist:

$\theta := \frac{d\alpha}{dx}$

Substituting in, the torque is found to be:

$\displaystyle T = \int_A rG(r\theta)dA = G\theta(\int_A r^{2} dA)$

Where: $\displaystyle dA = dydz = rdrd\theta$

The $2^{nd}$ polar moment of inertia for the beam is equal to the second term in the expression of the torque, $\displaystyle J = (\int_A r^{2} dA)$.

Fig 7. Comparison of two circular beams with equal area.

Thin-Walled Hollow Circular Cross-Section

Figure 7 shows two circular cross-sections, one solid and the other thin-walled, with the same area. The two cross-sections have different $2^{nd}$ polar moments of inertia.

For the thin-walled cross-section,

Inner radius: $\displaystyle r_i = a$

Outer radius: $\displaystyle r_o = b$

The $2^{nd}$ polar moment of inertia for the thin-walled cross-section can be determined through the expression:

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

Define $\displaystyle \bar{r}$ as the average radius:

$\displaystyle \bar{r} = \frac{a+b}{2}$

Define $\displaystyle t$ as the thickness:

$\displaystyle t = (a - b)$

Substitute into the $2^{nd}$ polar moment of inertia equation to get:

$\displaystyle J = \frac{1}{2}\pi(t)(2\bar{r})(b^{2} + a^{2}) = \pi(t)(\bar{r})(b^{2} + a^{2})$

Approximate and substitute:

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

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

$\displaystyle \Rightarrow J = \pi t \bar{r}(\bar{r}^{2} + \bar{r}^{2})= 2\pi t \bar{r}^{3}$

This result can be separated to allow for a proportional factor to be found between the $2^{nd}$ polar moment of inertia and the the average area, $\displaystyle \bar{A}$.

$\displaystyle \bar{A} = \pi\bar{r}^{2}$

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

In this form, it is clear that $\displaystyle J$ is proprtional to $\displaystyle \bar{A}^{\frac{3}{2}}$ by a factor of $\displaystyle (2\pi^{-\frac{1}{2}}t)$.

Torsion of a Uniform Non-Circular Bar

A uniform non-circular bar leads to warping of the cross-section. Warping is defined as the axial displacement along the x-axis (i.e. along bar length) of a point on the deformed rotated cross-section.

Fig 8. Torsion of Uniform Bars

Fig 9. 2-D view of Torsion on the Uniform Bar and the displacement vector PP'

Fig 10. Coordinate System show the Displacement in the y and z directions.

$U_y =$ y component of displacement vector $\vec{PP}'$

$U_z =$ z component of displacement vector $\vec{PP}'$

Important Definition : $\displaystyle \theta = \frac{\displaystyle \alpha}{x}$ = Rate of Twist

Referring to Figure 8, assuming that the cross-section at x=0 remains stationary, and if the rotation angle $\displaystyle \alpha$ is small, then the displacement components at P are given by the following:

$u = -r \displaystyle \alpha sin\displaystyle \beta = -\alpha y$ , which can be simplified to Equation (3.11):

(1) $U_y = - \displaystyle \theta x z$

Now the we have the displacement in the y direction, we can solve for the displacement $U_z$ which yields:

$\displaystyle U_z = + PP' cos \beta = + \alpha y_p$ , this can be simplified to Equation (3.12):

(2) $U_z = \displaystyle \theta x y$

The warping displacement in the x-axis direction can now be found as a function of the y and z directions.

(3) $U_x = \displaystyle \theta \displaystyle \psi (y, z)$

Based on these three equations come the Kinematic Assumptions.

Small Angle Approximation

---

Why is the $\vec{PP}'$ orthogonal to OP and P' is not on the circle, centered at O with a radius of OP?

Clearly from Figures 9 and 10, $\vec{OP}' > \vec{OP}$ and since the angle $\displaystyle \alpha$ is very small, the displacement can be though of as perpindicular to OP. Therefore, by the use of the small angle approximation for $\displaystyle \alpha$ , $OP = OP' = R$ . The length of displacement can also be reduced do the the small angle approximation. $U_z = R sin \displaystyle \alpha = R \alpha$ $U_y = R (1 - cos\displaystyle \alpha) = 0$

What is the Small Angle Approximation?

Fig 11: Use of the Small angle approximation for this case shows the value of the tangent of the angle is approximately just the value of that angle. Small angle approximation involves the simplification of the laws of trigonometry. The simplification is only approximately true for finite angles but becomes true as the limit of the angle approaches zero . Some simple examples include: $sin \displaystyle \theta = \theta$ $cos \displaystyle \theta = 1$

Fig 12: Airfoil with spars as shown.

Road Map for Torsional Analysis of Aircraft Wing:

---

A) Kinematic Assumption

B) Strain Displacement

C) Equilibrium Equation (stresses)

D) Prandtl Stress Function $\displaystyle \phi$

E) Strain Compatibility Equation

F) Equation for $\displaystyle \phi$

G) Boundary Conditions for $\displaystyle \phi$

H) $T = 2 \iint_A \displaystyle \phi \,dA$

$T = G J \displaystyle \theta$

$J =\frac{-4}{\nabla^2 \phi}\iint_A \displaystyle \phi \,dA$

I) Thin-Walled Cross-Section , $T = 2 q A$ (3.48) This equation was also derived in Problem 1.1.

J) Twist Angle $\displaystyle \theta$ Method 1.

$\displaystyle \theta = \frac{1}{2 G A} \oint \frac{q}{t}\,ds$

Proving the Area of a Triangle is $\frac{b\ h}{2}$

In the lecture on the 14th lecture on the 26th of September, 2008 the professor mentioned how understanding the fundamentals is important and assigned the homework problem of proving the area of a triangle.

To begin this proof, it will be assumed that the area of a parallelogram $A_p$ is understood and given to be

Fig 13. Parallelogram split in half into two triangles.

$A_p = b\ h$

To begin, if it is stated that

$\bar {AB} \parallel \bar {CD}$

$\bar {BC} \parallel \bar {AD}$

It is understood that the given shape is a parallelogram with an area given above. Since it is a parallelogram, it is understood that

$\bar {AB} \cong \bar {CD}$

$\bar {BC} \cong \bar {AD}$

And since these sides are congruent, by using the side angle side theorem of geometry (SAS) it can be understood that

$\angle ADB \cong \angle CBD$

$\angle ABD \cong \angle BCD$

Hence by the angle side angle theorem of geometry (ASA) it is understood that

$\triangle ABD \cong \triangle BCD$

Since the parallelogram is made up of these two congruent triangles, it is understood that

$A_{t} = \frac{b\ h}{2}$

It can be noted that any triangle given can be extended into one or two parallelograms, dependent on the triangle, and superimposed to have the result that was proven above.

QED

Comparison of the Polar Moment of Inertia of a Solid Circle and a Thin walled circle

The polar moment of inertia is used to quantify the ability of an object to resist torsion. It can be calculated by the definition

$J_z = \iint_A \rho^2\ dA$

Where $\rho$ is the radial distance from the z axis to the element $dA$. $dA$ is an infinitesimal area element. Converting $dA$ to polar coordinates which was done in ***homework 2*** the term $dA$ expands to

$dA = \rho\ d\rho\ d\theta$

Solid Circle Proof

When substituting the above expressions, the polar moment of inertia of a solid circle becomes:

$J_z = \int_{0}^{2\pi} \int_{0}^{R} \rho^3\ d\rho\ d\theta$

and after integrating results with

$J_z = \frac{\pi\ R^4}{2}$

Where $R$ is the radius of the circle.

Thin walled Circle Proof

As is possible with the second area moment of inertia, its possible to add and subtract the polar moment of inertia from others which result in the desired shape. Since the polar moment of inertia of a solid cross section has already been proven, it would be prudent to use this to our advantage by subtracting the polar moment of inertia of the inner "hollow" section from the polar moment of inertia of the outer shelled section.

$J_{z,out} = \frac{\pi\ R_o^4}{2}$

$J_{z,in} = \frac{\pi\ R_i^4}{2}$

$J_z = J_{z,out} - J_{z,in} = \frac{\pi}{2}\ (R_o^4 - R_i^4)$

Fig 14. Circle showing that $\rho$ is the distance to the differential element dA.

Now if we make the assumption that $t \ll R_o$. Noting that

$R_o \approx R_i + t$

It makes it possible to define a new variable $\bar r$ which is defined as

$\bar r := \frac{R_o + R_i}{2} = \frac{2\ R_o - t}{2} = \frac{2\ R_i + t}{2}$

Now assuming that $t \ll R_o$ it can be approximated that

$\bar r^2 = (R_o - \frac{t}{2})^2 = R_o^2 - R_o\ t + \frac{t^2}{4} \approx R_o$

$\bar r^2 = (R_i + \frac{t}{2})^2 = R_i^2 + R_i\ t + \frac{t^2}{4} \approx R_i$

To the previous equation of the polar moment of inertia

$J_z = \frac{\pi}{2}\ (R_o^4 - R_i^4) = \frac{\pi}{2}\ ((R_o^2)^2 - (R_i^2)^2)$

can be expanded to

$J_z = \frac{\pi}{2}\ (R_o - R_i)\ (R_o + R_i)\ (R_o^2 + R_i^2)$

noting that with the approximations it can be given that

$(R_o - R_i) = t$

$(R_o + R_i) = 2\ \bar r$

$(R_o^2 + R_i^2) = 2\ \bar r^2$

and substituting these approximations

$J_z = \frac{\pi}{2}\ (t)\ (2\bar r)\ (2\bar r^2) = 2\ \pi\ t\ \bar r^3$

Comparing Solid versus thin walled

As shown above, the polar moment of inertia of a solid circle is

$J_{z,solid} = \frac{\pi\ R^4}{2}$

and for a thin walled is

$J_{z,thin} = 2\ \pi\ t\ \bar r^3$

When analyzing these equations it can be shown that most of the torsional resistance is a result of the outer material rather than the inner material. This means that the material at the center of a solid circle is essentially wasted in comparison to the outer material. If the same amount of material were used, the thin walled geometry gives a lot more bang for its buck.

In applications such as aerospace structures where weight is an extremely large factor, less material is always better. This would result in the use of thin walled circular cross-section for use in structural integrity throughout the aircraft.

Problem 1.1 Applied to General Single Cell Section

Commonly in aerospace applications it may be necessary to determine the characteristics of a hollow, thin-walled cross-section similar in shape to that of an airfoil. In this case, the cross-section can be subdivided into smaller parts, namely, one triangular and one semi-circular cross-section. Each respective $\bar {A}$ can be calculated and summed in order to determine the total $\bar {A}$ of the section as follows:

$\bar {A_{circ}} = \frac{\pi}{2}(\frac{b}{2})^2$ and $\bar {A_{tri}} = \frac{b a}{2}$. Therefore,

$\bar {A_{total}} = \frac{\pi}{2}(\frac{b}{2})^2 + \frac{b a}{2} = \frac{\pi}{2}(\frac{2}{2})^2 + \frac{2(4)}{2} = \frac{\pi}{2} + 4 = 5.5708 m^2$

Fig 15 - Problem 1.1 applied to an airfoil-like single cell section

The torque generated by a single cell is: $T_i = 2q_i \bar {A_i}$

while the total torque is: $T = \sum_{i=1}^{n_c} T_i = \sum_{i=1}^{n_c} 2q_i \bar{A_i}$ where $n_c$ is the total number of cells.

Also, if the value of shear flow is desired, the above equations could be rearranged to yield:

$shear flow = q = \frac{T}{2 \bar {A}}$

Additionally, the twist angle could be calculated as:

$twist angle = \theta = \frac{1}{2 G \bar {A}} \sum_{j=1}^{n_c} \frac{q_j l_j}{t_j}$

where G is the shear modulus and $n_c$ is the number of cells.

Quadrature

The concept of quadrature (also known as squaring the circle) has its roots in ancient geometry. The problem consists of constructing a square from a circle of equal area using only a compass and straightedge and accomplishing this task in only a finite number of steps. The problem dates back to the 19th century BCE and has been attempted by numerous widely-known mathematicians. Although many approximations of the solution exist, it has been proven that a solution to the problem is impossible. It was proven in the early 19th century CE that if $\pi$ were a transcendental, then a solution to the problem does not exist. In 1882, Ferdinand von Lindemann proved in his Lindemann-Weierstrass Theorem that $\pi$ is, in fact, a transcendental, therefore, proving that squaring the circle is impossible.

Fig 16 - The circle and square have equal areas. The concept of Squaring the Circle entails constructing the square from the circle shown.

MATLAB Question

The task for the MATLAB portion of the HW report involves plotting an arbitrary NACA airfoil, finding the average area of the airfoil, and locating the airfoil's area centroid.

Full MATLAB Code

function abar(bool,c,m,p,t,ns,P_o) %Matlab Function %************************************** %filename: abar.m % %Purpose: % Determine the average area of an arbitrary NACA 4-digit series airfoil % and also calculate the centroid of the airfoil section. % %Author: % Team VQCrew (eas4200c.f08.vqcrew) % EAS4200C (Aerospace Structures) Spring 2008 % % Modified on: Wed, 08 Oct 2008 % Created on: Mon, 06 Oct 2008 % %Remarks: % % Command line syntax: % abar(bool,c,m,p,t,ns,P_o) % % Variables: % bool: Choose either a circular (0) or standard (1) airfoil (boolean) % c: Desired chord length (real) % m: Maximum camber in percentage chord (real) % p: Position of maximum camber in tenths chord (real) % t: Maximum thickness in percentage chord (real) % ns: Number of segments for contour subdivision (real) % P_o: Reference position for area calculation (column vector - [y;z]) % % Notes: % All of the variables listed above must be provided when calling the % function for proper program operation. % P_o must be input as a column vector representing the coordinates of the % reference point. % % Three subfunctions called in the program: % [yu,zu,yl,zl]= airfoil(c,m,p,t,ns) - generates standard airfoil contour % [yu,zu,yl,zl] = cirairfoil(c,ns) - generates circular airfoil contour % [cent]=centroid(yu,zu,yl,zl,n) - calculates centroid of area %/<--function start-->/% % select circular or regular airfoil contour if bool == 1 [yu,zu,yl,zl] = airfoil(c,m,p,t,ns); else [yu,zu,yl,zl] = cirairfoil(c,ns); end n = length(yu); % iterate counter-clockwise over the upper surface (trailing edge to % leading edge) using vector cross-product method for area calculation for j=n:-1:2 Pu = [yu(j);zu(j)]; ru = Pu - P_o; Qu = [yu(j-1);zu(j-1)]; ru = [ru;0]; lu = Qu - Pu; lu = [lu;0]; Bu = .5*cross(ru,lu); Au(j) = Bu(3); end % iterate counter-clockwise along the lower surface (leading edge to % trailing edge) using vector cross-product method for area calculation for i=1:n-1 Pl = [yl(i);zl(i)]; Ql = [yl(i+1);zl(i+1)]; rl = Pl - P_o; rl = [rl;0]; ll = Ql - Pl; ll = [ll;0]; Bl = .5*cross(rl,ll); Al(i) = Bl(3); end % calculate airfoil cross-section area A_bar A1 = sum(Au); A2 = sum(Al); A_bar = A1 + A2; % calculate the centroid of the given airfoil [cent] = centroid(yu,zu,yl,zl,n); % print relevant data to command window fprintf('\n') fprintf('The average area of the airfoil is %1.3f square chord units\n',A_bar) fprintf('\n') fprintf('The centroid is located at (%1.3f,%1.3f) chord units\n',cent(1),cent(2)) fprintf('\n') % plot airfoil contour and centroid cross-hair plot(yu,zu,'-b',yl,zl,'-b') hold on plot(cent(1),cent(2),'k+') axis([0 c -.5*c .5*c]) end %***************************************** function [yu,zu,yl,zl]=airfoil(c,m,p,t,ns) % generate subdivided y-axis points y = linspace(0,1,ns); n = length(y); % differentiate between cambered and symmetric airfoils if p==0 || m==0 z = zeros(1,n); dz = zeros(1,n); t = t/100; % convert thickness to percentage chord else % convert NACA 4 digit code into percentage chord/tenths chord % note, chord length will be adjusted for later m = m/100; p = p/10; t = t/100; % calculate mean camber line contour and slope, pre-max camber i=1; while 1 if y(i) >=p, break, end z(i) = (m./p.^2).*(2.*p.*y(i) - y(i).^2); dz(i) = (m./p.^2).*(2.*p - 2.*y(i)); i = i+1; end % calculate mean camber line contour and slope, post-max camber for j = i:n z(j) = (m./((1-p)^2)).*((1-2.*p)+2.*p.*y(j) - y(j).^2); dz(j) = (m./((1-p).^2))*(2.*p - 2.*y(j)); end end % calculate thickness distribution ta = abs(5.*t.*(.2969.*sqrt(y)-.1260.*y-.3516.*y.^2 + .2843.*y.^3 - .1015*y.^4)); % calculate mean camber line tangents theta = atan(dz); % calculate airfoil contour coordinates, with adjustment for given chord % length yu = (y - ta.*sin(theta))*c; zu = (z + ta.*cos(theta))*c; yl = (y + ta.*sin(theta))*c; zl = (z - ta.*cos(theta))*c; end %************************************** function [yu,zu,yl,zl]=cirairfoil(c,ns) %generate circular airfoil contour coordinates y = linspace(0,1,ns); ys = y-.5; z = abs(sqrt(.25-ys.^2)); yu = y*c; yl = y*c; zu = z*c; zl = -z*c; end %************************************** function [cent]=centroid(yu,zu,yl,zl,n) asum = 0; cen_y = 0; cen_z = 0; % calculate centroid of airfoil using trapezoidal Riemann sum for i=1:n-1 h_u = yu(i+1)-yu(i); h_l = yl(i+1)-yl(i); bs_u = zu(i+1)+zu(i); bs_l = abs(zl(i+1)+zl(i)); a_u = (h_u*bs_u)/2; a_l = (h_l*bs_l)/2; asum = asum + a_u + a_l; cenu_y = h_u*(2*zu(i+1)+zu(i))/(3*(zu(i+1)+zu(i))) + yu(i); cenl_y = h_l*(2*zl(i+1)+zl(i))/(3*(zl(i+1)+zl(i))) + yl(i); cenu_z = (zu(i)^2 + zu(i)*zu(i+1) + zu(i+1)^2)/(3*(zu(i+1)+zu(i))); cenl_z = (zl(i)^2 + zl(i)*zl(i+1) + zl(i+1)^2)/(3*(zl(i+1)+zl(i))); cen_y = cen_y + (cenu_y*a_u) + (cenl_y*a_l); cen_z = cen_z + (cenu_z*a_u) + (cenl_z*a_l); end cent(1) = cen_y/asum; cent(2) = cen_z/asum; end

Using the following input (in accordance with the syntax presented in the code header), the code generates Figure 17, a graph of the NACA 2415 airfoil with a cross-hair at the centroid. The output references "chord-units", as the code was purposely generated to handle any unit system for chord length. For this particular case, the chord length is set at .5 m.

abar(1,.5,2,4,15,1000,[0;0])

The average area of the airfoil is 0.025 square chord units

The centroid is located at (0.210,0.008) chord units

Fig 17 - MATLAB-plotted NACA 2415 airfoil

Contributing Team Members

The following students contributed to this report:

John Saxon Eas4200c.f08.vqcrew.c 20:40, 8 October 2008 (UTC)
Javier Stober Eas4200c.f08.vqcrew.b 14:25, 8 October 2008 (UTC)
Darin Toscano Eas4200c.f08.vqcrew.d 17:58, 8 October 2008 (UTC)
Kevin Klauk Eas4200c.f08.vqcrew.E 19:34, 8 October 2008 (UTC)
Steven Hepsworth Eas4200c.f08.vqcrew.f 17:57, 8 October 2008 (UTC)
Adam Edstrand Eas4200c.f08.vqcrew.a 20:20, 8 October 2008 (UTC)