User:Eas4200c.fo8.aero.lee/HW6

 See my comments below. Please don't remove these comment boxes; just add your comment to these comment boxes in case you fix the problems.

Warning: User Eas4200c.f08.lee was not registered; please register to create this user account, since otherwise the page may be deleted.

Also, don't put the figures on the right of the collapsible tables, since these figures would block the link to expand/collapse these tables; see your figures to the right of your matlab codes in collapsible tables below. Instead, you may want to put these figures to the left of the collapsible tables. Another way is to shorten the width of the collapsible tables.

Eml4500.f08 13:27, 23 November 2008 (UTC)

=Plate Buckling= A plate under uniform stress in one dimension (see Figure 1) undergoes deformation according to the following equation:
 * $$u_z(x,y)= \sum_{m=1}^{ \infty }{}\sum_{n=1}^{\infty}\psi_{mn} (x,y) $$ where $$\psi$$ is expressed as
 * $$\psi_z(x,y)=csin(\frac{m\pi x}{a})sin(\frac{n\pi x}{b})$$

and $$u$$ is the out of plane deformation. (Can be thought of as deformation in z direction.) Please note this is for a simply supported plate.

The resulting shape when looking at the plate normal to the xy plane is called the buckling mode shape.

The critical loading, $$P_x$$ is given by the following equation:
 * $$(P_x)=k_c\frac{\pi ^{2}D}{b}$$


 * $$k_c(m,\frac{a}{b})=(\frac{mb}{a}+\frac{a}{mb})^2$$

and


 * $$D=\frac{Eh^3}{12(1-\upsilon ^2)}$$

Where h is the thickness of the plate in the z direction.

The critical buckling stress is:
 * $$(\sigma_{xx})_{cr}=\frac{(P_x)_{cr}}{bh}=k_c \frac{\pi^2D}{b^2h}$$

Components related to the minimum critical value are denoted with a star, e.g. $$\sigma^\star_{cr}$$ The minimum critical buckling stress is dependent on $$k^\star_c$$, which in turn is dependent on $$m^\star$$ by the following equation:
 * $$k_c^\star(a/b)=\min_{m = 1, 2, \cdots}\ k_c(m, a/b)=\min_{m = 1, 2, \cdots}\left( \frac {m b}{a} +\frac{a}{m b}\right)^2 =

\left( \frac{m^\star b}{a}+\frac{a}{m^\star b}\right)^2$$ where
 * $$ m^\star(a/b)=\underset{m=1,2, \cdots}{\rm argmin}\ k_c(m, a/b) =\underset{m=1,2, \cdots}{\rm argmin}\left(\frac{m b}{a}

+ \frac {a}{m b}\right)^2$$ The transition from m to m+1 waves occurs when:
 * $$\frac{a}{b}=\sqrt{m(m+1)}$$

For a square plate with all 4 edges clamped, the buckling mode shape approximates to:
 * $$ \psi(x,y)=\frac{c}{4}\left(1 - \cos \frac{2 \pi x}{a}\right)\left(1 - \cos \frac{2 \pi y}{b} \right)$$

This gives a “bubble” like shape, as seen in the matlab figure below. It makes logical sense since the slope of the buckled shape must approach zero (dz/dx and dz/dy) at all clamped edges.

Shear Load, Simply Supported
Load is denoted as $$N_{xy}$$, which is force per unit length, or stress times h, or plate thickness.

Stored strain energy is equated to the work done by shear load U=W.

The Nxy term is in the equation for strain energy, and the critical Nxy is expressed as:
 * $$\displaystyle

(N_{xy})_{cr} =-\frac{ab D}{32}\frac{\displaystyle\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}C_{mn}^2\left(\frac{m^2 \pi^2}{a^2}+       \frac{n^2 \pi^2}{b^2}\right)^2}{\displaystyle\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}\sum_{p=1}^{\infty}\sum_{q=1}^{\infty} C_{mn}C_{pq}\frac{mnpq}{(m^2 - p^2)(n^2 - q^2)}} $$ After further analysis and derivation, crit. buckling stress is:
 * $$(\sigma_{xy})^{\star,5}_{cr}=\pm k^{\star, 5}_c \frac{\pi^2D}{b^2h}, k^{\star,5}_c = \frac{\pi^2}{32}\frac{1}{\lambda\vartheta }$$

=A note on dimensional analysis=

A bracket notation (ex. [m]) is used to indicate the dimensions associated with a variable. Here are some examples:


 * $$[f]=\frac{F}{L}$$

this means the dimensions of "f" are force per length.

Since


 * $$[A]=L^2$$

then


 * $$\left[\frac{f}{A} \right]= \frac{F}{L^{3}}$$

meaning the dimensions of "f" per unit area are force per volume, or length cubed. Note, L is only a unit of measure, not a specific dimension.

The following describes a body force:


 * $$\frac{d\sigma }{dx}+\frac{f(z)}{A}=0\Rightarrow [\sigma ]=\frac{F}{L^2}\Rightarrow \left[\frac{d \sigma}{dx} \right]= \frac{[d \sigma ]}{[dx]}=\frac{F/L^2}{L} = \frac{F}{L^3}$$

To show that axial strain is nondimensional:


 * $$[\epsilon ] = \frac{du}{dx} = \frac{L}{L} = 1...nondimensional...$$

Shear strain is also nondimensional:


 * $$[\gamma ]=\frac{[\epsilon _{yy}]}{[\epsilon _{xx}]}=1...nondimensional...$$

Also note
 * $$\left[\frac{d \sigma_{ij}}{dx_i} \right]=\frac{F}{L^3}=\frac{Force}{Volume}$$

Intro
We introduce a function, $$\phi$$, which plays the role of a potential function. We call this the Prandtl stress function. $$\sigma_{yx}$$ and $$\sigma_{zx}$$ are components of the gradient of $$\phi$$ with respect to y and z.

In wikipedia, a gradient is defined as is a vector field which points in the direction of the greatest rate of increase of the scalar field, and whose magnitude is the greatest rate of change. In notation, it is defined as
 * $$\bigtriangledown f(x,y,z)= \frac{\partial f}{\partial x}i + \frac{\partial f}{\partial y}j + \frac{\partial f}{\partial z}k$$

where f is the potential.

For our case,
 * $$\sigma_{yx}=\frac{\partial \phi}{\partial z}$$ and
 * $$\sigma_{zx}=-\frac{\partial \phi}{\partial y}$$

We can plug these values into the following equation:
 * $$\frac{\partial \sigma_{yx}}{\partial y}+ \frac{\partial \sigma_{zx}}{\partial z}=0$$

to obtain
 * $$\frac{\partial }{\partial y}\frac{\partial \phi}{\partial z}+\frac{\partial }{\partial z} \frac{(-\partial \phi)}{\partial y}=\frac{\partial^2 \phi}{\partial y \partial z}- \frac{\partial^2 \phi}{\partial z \partial y}=0$$

Phi is continuous and smooth, making the 2nd mixed derivative interchangeable, i.e.
 * $$\frac{\partial^2 \phi}{\partial y \partial z}= \frac{\partial^2 \phi}{\partial z \partial y}$$

These equations can be compared with Mechanics of Aircraft Structures(Sun) Ch. 3 eqs 3.15-3.19. Using the author's notation:
 * $$\frac{\partial \tau_{yz}}{\partial x}-\frac{\partial \tau_{xz}}{\partial y}=2G\theta$$

and
 * $$\frac{\partial^2 \phi}{\partial x^2}+ \frac{\partial^2 \phi}{\partial y^2}=-2G\theta$$

The remainder of section 3.2 brings us to our roadmap equation for J (Roadmap eq 7) and notes that once the Prandtl stress function is solved, GJ, the torsional rigidity (the bar's ability to resist deformation due to torsion) is also found. Torsional analysis of structures fall into two cases - that of structures with a thin walled cross section (also referred to as closed cross section) and those of structures with a solid cross section.

Case 1: Thin Walled Cross Section
A prime example of this case is that of the NACA airfoils we have been studying throughout the semester. For thin walled cross section the boundary conditions for Prandtl stress function, $$ \phi $$, are that $$\phi$$ is a constant on inner and outer surfaces of the cross section. The figure below depicts the inner and outer surfaces, $$ S_1 $$ and $$ S_0 $$, respectively.



The boundary conditions for this case are $$\phi=C_0$$ on $$S_0$$ and $$\phi=C_1$$ on $$S_1$$.

Case 2: Solid Cross Section
This case consists of any bar or beam that has a solid cross section such as the one depicted below.



The constant value of $$\phi$$ at the boundary solid cross sections is arbitrary and so it can be set to zero on that surface. Thus the boundary in this case is $$\phi = 0$$ on $$S_0$$

Uniform Bar With Solid Circular Cross Section
As an example, we will analyze the following bar with a circular cross section as shown below.



For this case, $$\phi$$ can be defined as


 * $$\phi(y,z) = C\left(\frac{y^2}{a^2}+\frac{z^2}{b^2}-1\right)$$

It might also be noted that for the case of a circular cross section a=b. As mentioned before for the case of a solid cross section $$\phi=0$$ on $$S_0$$. Recalling the relation:


 * $$\nabla^2\phi=-2G\Theta$$

an expression for C can be obtained by solving for the laplacian of $$\phi$$.


 * $$\nabla^2\phi = \frac{4C}{a^2} = -2G\Theta$$
 * $$\Rightarrow C = -\frac{1}{2}a^2G\Theta$$

To analyze the Prandtl stress function, $$ \phi $$, as a function of the radius of a circular cross section of diameter $$ 2\alpha $$, consider the definition of the radius as


 * $$ r = \left(y^2 + z^2 \right)^{1/2} $$

Then, a plot of $$ \phi(r) $$ vs. r would yield the figure below



Now, recall from previous lectures the following equations


 * $$ T = 2\int_{A}^{}{\phi dA} = 2C\left(\frac{J}{a^2}-A\right) $$


 * $$ J = \int_{A}^{}{r^2 dA} = \frac{1}{2} \pi a^4 $$


 * $$ \tau = \frac{Tr}{J} $$

Combine these equations with the following


 * $$ \sigma_{yx} = \frac{\partial \phi}{\partial z} = -G \theta z $$


 * $$ \sigma_{zx} = -\frac{\partial \phi}{\partial y} = G \theta y $$

to obtain


 * $$ u_x\left(y,z\right) = 0 $$

This implies that there is no axial displacement at any point along a bar with uniform, solid, circular cross-section. This was assumed to be true during ad hoc derivations of twist; however, we see now that the Theory of Elasticity proves it to be so. A rigorous proof is presented below.

 Error: See my comments in Team VQCrew (with annotation). Eml4500.f08 14:04, 23 November 2008 (UTC)

= Traction Force Equation: Derivation =

Geometric Manipulation Using an Infinitesimal Element
Goal: Using the traction force equation to derive the equation:

$$\displaystyle[t_i]_{3x1}=[\sigma_{ij}]_{3x3}[n_i]_{3x1} $$

Consider the following figure below of an infinitesimal element to be able to be able to derive geometrically, the traction force from the stress tensor and normal vector at the surface of the element.



From the figure the vector $$\boldsymbol\vec{n}$$ has the vector equation: $$\displaystyle \boldsymbol \vec{n}=n_y \boldsymbol \vec{j}+n_z \boldsymbol \vec{k}$$

From this equation, $$\boldsymbol\vec{n}$$ is the unit normal vector and is given as $$\displaystyle ||\boldsymbol \vec{n}||=1$$

Looking at the figure, the elemental geometry helps in the establishing the following relationship between dy and dz:

$$\displaystyle dy=ds*sin(\theta)$$ $$\displaystyle dz=ds*cos(\theta)$$

Additionally, the normal vector equation can be used to rewrite the vector in terms of angle $$\theta$$ and conclude also that:

$$\displaystyle n_y=cos(\theta)$$ $$\displaystyle n_z=sin(\theta)$$

Looking at Figure and continuing through this derivation, the next logical choice would be to take the sum the forces in the y and z directions using the equilibrium equations. It is important to note that the depth into the page is taken as the x-direction and is neglected for this case since the depth is considered to be unity, or 1.

As given in class, the derivation for the equilibrium in the y direction of the sum of forces is given as:

$$\sum {F}_y=0=-\sigma_{yy}*(dz*1)-\sigma_{yz}*(dy*1) +t_y*(d_s*1)$$

Substituting values for dy and dz into the equilibrium equation above, the following equation is produced:

$$\sum {F}_y=0=-\sigma_{yy}ds n_y-\sigma_{yz}ds n_z +t_yd_s$$

Further simplification of the equilibrium equation above and solving for $$\displaystyle t_y : $$

$$\displaystyle t_y= \sigma_{yy}n_y+\sigma_{yz}n_z(1)$$

As part of HW 6, the equilibrium equation for the sum of forces in the z-direction are given below:

Using equations (1) and (2) for tractional force components, we can arrange the equations in the form of a matrix, reflected below to produce equation 3:

$$\begin{Bmatrix} t_y\\ t_z\\ \end{Bmatrix} = \begin{bmatrix} \sigma_{yy} & \sigma_{yz}\\ \sigma_{zy} & \sigma_{zz}\\ \end{bmatrix} \begin{Bmatrix} n_y\\ n_z\\ \end{Bmatrix}(3) $$ Similarly for a general 3D case, a similar matrix may be constructed using indicial notation, shown below:

$$\begin{Bmatrix} t_1\\ t_2\\ t_3\\ \end{Bmatrix} = \begin{bmatrix} \sigma_{11} & \sigma_{12}& \sigma_{13}\\ \sigma_{21} & \sigma_{22}& \sigma_{23}\\ \sigma_{31} & \sigma_{32}& \sigma_{33}\\ \end{bmatrix} \begin{Bmatrix} n_1\\ n_2\\ n_3\\ \end{Bmatrix}(4) $$

An advantage of using this indicial notation allows us to rewrite the entire matrix equation in the form of a single equation that is more compact. This compact notation is shown below:

$$t_i=\sum_{j=1}^3{\sigma_{ij}n_j},i=1,2,3$$

In this equation looking at the stress value, where the subscript i, and j designates the row, and column of the matrix respectively.

Generalizing this formula, produces our end goal and solves for the desired equation: $$\displaystyle [t_i]_{3x1}=[\sigma_{ij}]_{3x3}[n_i]_{3x1} $$

Roadmap Incorporation
Looking at the roadmap on lecture page 16-2, we recall the following three equations:

$$\displaystyle T=2\int\int_{A}^{}{}\phi dA(5)$$ $$\displaystyle T = GJ\theta~(6) $$ $$\displaystyle J=\frac{-4}{\nabla^{2}\phi }\int \int _{A}\phi dA(7)$$

These three equations can be combined, and the equation (8) shown below can be derived by solving the three equations shown above for the Lagrangian of phi.

$$\frac{\partial^2{\phi}}{\partial{y}^2}+\frac{\partial^2{\phi}}{\partial{z}^2}=-2G\theta(8)$$

This becomes important because $$\phi$$ is taken to be constant on the lateral surface of the bar.

= Flexural Shear in Thin-walled Sections =

Motivation
It is known that for a generic loading condition of a beam, there exists not only a bending moment but also a transverse shear load. Consider the beam in the figure below. Some generic loading condition is applied to the straight beam on the left. A cut can be made in the beam at a distance x from the wall. The bending moment and shear loads are drawn for each side of the exposed section, where M(x) represents the moment as a function of distance x and V(x) represents the shear load as a function of distance x.



= Unsymmetrical Thin-walled Cross Section =

General Equation for unsymmetrical cross sections $$\displaystyle \int_{A_s} \frac{\sigma_{xx}}{dx} dA = -q_s $$ Eq. 5.1

Symmetric Cross Section (about y- axis)
$$ \displaystyle \sigma_{xx} = \frac{M_y z}{I_y}$$

$$ \displaystyle q(s) = \frac{-V_z Q_y}{I_y}$$ Eq. 5.2

$$ \displaystyle Q_y = \int_{A_s} z dA = z_c A_s$$

Unsymmetric Cross Section
$$ \displaystyle \sigma_{xx} = (k_y M_z - k_{yz} M_y)y + (k_z M_y - k_{yz} M_z)z$$ → Eq. 5.3

This is true for the following k values from p.26-2 of notes: $$ \displaystyle k_y = \frac{I_y}{D}$$ $$k_{yz} = \frac{I_{yz}}{D}$$ $$k_z = \frac{I_z}{D}$$

The expression for $$\displaystyle \sigma_{xx}$$ can also be put into matrix form:

$$\displaystyle \sigma_{xx} = \begin{bmatrix} z & y\\ \end{bmatrix} \begin{bmatrix} k_z & -k_{yz}\\ -k_{yz} & k_y\\ \end{bmatrix} \begin{Bmatrix} M_y\\ M_z\\ \end{Bmatrix} $$

$$ \displaystyle q(s) = -(k_y V_y - k_{yz} V_z) Q_z - (k_z V_z - k_{yz} V_y) Q_y$$ → Eq. 5.5

$$ \displaystyle => Q_z = \int_{A_s} y dA$$ and $$ \displaystyle Q_y = \int_{A_s} z dA $$

= Matlab Analysis =



The main purpose of this weeks homework problem is to examine the bending stress in the skin of the NACA airfoil. In order to determine this, we will use exact approximations as well as approximations for simply supported and clamped boundary conditions since the airfoil's spars and stringers cause the skin to fall somewhere in between the two.

The first part involves using the previous homework's inertia matrix, centroid, and applied moments to determine the normal bending stress in the skin between the two spars as a function of the y-coordinate using the equations:


 * $$\sigma_{xx}=E\epsilon_{xx}=\frac{I_yM_z-I_{yz}M_y}{I_yI_z-(I_{yz})^2}y+\frac{I_zM_y-I_{yz}M_z}{I_yI_z-(I_{yz})^2}z$$



\begin{bmatrix} I_{22} \\ I_{33} \\ I_{23} \\ \end{bmatrix}=

\begin{bmatrix} 1.298\ x\ 10^{-5}\ m^4 \\ 6.346\ x\ 10^{-7}\ m^4 \\ -9.494\ x\ 10^{-8}\ m^4 \\ \end{bmatrix}$$


 * $$Centroid(y,z)=(0.210\ m,\ 0.210\ m)$$

 The figure does not look right: The normal stress $$\displaystyle \sigma_{xx}$$ is expected to rise and fall with the coordinate $$\displaystyle z$$ of the airfoil skin, i.e., the function would look curved similar to the curve of the airfoil skin, instead of a straight line. See, e.g., Team VQCrew (with annotations) Eml4500.f08 11:36, 26 November 2008 (UTC)



From this we find that the average stresses are $$\sigma_{xx,FB}=26.5\ MPa$$ and $$\sigma_{xx,EH}=32.1\ MPa$$.

Next we wish to explore what shapes the skin would undergo due to these bending stress. For an aspect ratio of 1.5, we can see in Figures XXX and XXX that depending on m and n, different modes of buckling occur. The variables m and n determine the number of half wave lengths in the x and y direction respectively, which is readily apparent in the figures.

We are interested in finding the period T of this function that determines the shape of the buckled plate to prove the m is indeed the number of half waves. This involves finding the period of $$ \sin ( \frac{m \pi x}{a} )$$  by considering   $$\sin ( \frac{m \pi (x + T)}{a} ) = \sin ( \frac{m \pi x}{a} )$$. We know that:


 * $$\sin( \theta ) = \sin( \theta + 2\pi k)=\sin(\frac{m \pi x}{a}+\frac{m \pi T}{a})$$


 * $$Period=T=\frac{2\pi}{\frac{m\pi}{a}}=\frac{2a}{m}$$

For constant a=1, we see that the period is 2/m. For m=1, we see the period is 2 and therefore there is one half of the wave at a=1. For m=2, the period is 1 and two half waves at a=1, so on and so forth.



 I moved this figure to the left, so the figure does not block the link to show the collapsible table containing your matlab code, and added a few blank lines, so this collapsible table does not get pushed too far to the right. Eml4500.f08 11:36, 26 November 2008 (UTC)

In the determining of critical stresses, a variable of large importance is k_c. This number is a function purely of m, the number of half wave lengths, and the aspect ratio of the plate. The figure to the right shows how k_c varies with these variables. This factor is used in critical stress calculations for simply supported plates.



 I moved this figure to the left, so the figure does not block the link to show the collapsible table containing your matlab code, and added a few blank lines, so this collapsible table does not get pushed too far to the right. Eml4500.f08 11:36, 26 November 2008 (UTC)

As we discussed earlier, there are stringers and spars connect to the skin so it is not just simply supported. The other approximation for bending of a flat plate is clamped boundary condition. Here the two ends are not allowed to rotate, making the slope of plate zero at the boundary. The figure to the right shows how this sort of boundary condition affects the buckling mode of the plate.



The final section is interested in modeling the skin section in between the two spars as flat plates and determining the critical stresses in this plate. We will compare the value reached earlier for average stress to each of the two boundary conditions. In this figure, we vary the aspect ratio of plate while holding the width constant to that of the actual skin section. The two cases shown are for the clamped and simply supported boundaries mentioned earlier. We can see that at small aspect ratios, our average stress for the top plate, which the figure was based upon, falls in between the two curves. At higher aspect ratio's however (which is what we would expect for an aircraft wing), the stress is higher the than the critical stress, implying that the skin does not work well to resist bending, which is what we expected.

= Contributing team members =


 * Gonzalo Barcia--Eas4200c.f08.aero.barcia 03:05, 21 November 2008 (UTC)
 * William Kurth--Eas4200c.f08.aero.kurth 14:56, 21 November 2008 (UTC)
 * Melisa Gaar--Eas4200c.f08.aero.gaar 19:21, 21 November 2008 (UTC)
 * Thomas McGilvray--Eas4200c.f08.aero.mcgilvray 19:38, 21 November 2008 (UTC)
 * Ray Strods--Eas4300c.f08.aero.strods 19:19, 21 November 2008 (UTC)
 * Jared Lee--Eas4200c.f08.aero.lee 19:43, 21 November 2008 (UTC)
 * Oliver Oyama--Eas4300c.f08.aero.oyama 19:54, 21 November 2008 (UTC)