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



Team member eas4200c.f08.vqcrew.b (Stober) forgot to sign the page when it was originally due in early November. The situation has been rectified. Eas4200c.f08.vqcrew.c 19:50, 10 December 2008 (UTC)

Linear Stress-Strain Relations
The expression giving the strain induced by a stress can be rewitten to allow for the deriviation of a second expression relating the stress induced for a given strain.
 * $$ \begin{bmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{33} \\ \varepsilon_{23} \\ \varepsilon_{13} \\ \varepsilon_{12} \end{bmatrix} = \begin{bmatrix} \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{bmatrix} \begin{bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{33} \\ \sigma_{23} \\ \sigma_{13} \\ \sigma_{12} \end{bmatrix} $$

Becomes,


 * $$\left\{\varepsilon_{ij} \right\}_{6x1} = \begin{bmatrix}A_{3x3} & 0_{3x3} \\ 0_{3x3} & B_{3x3} \end{bmatrix}\left\{\sigma_{ij} \right\}_{6x1}$$

Where $$A_{3x3} = \begin{bmatrix} \frac{1}{E} & \frac{-\nu}{E} & \frac{-\nu}{E} \\ \frac{-\nu}{E} & \frac{1}{E} & \frac{-\nu}{E} \\ \frac{-\nu}{E} & \frac{-\nu}{E} & \frac{1}{E}\end{bmatrix}$$ and $$B_{3x3} = \begin{bmatrix}\frac{1}{2G} & 0 & 0 \\ 0 & \frac{1}{2G} & 0 \\ 0 & 0 & \frac{1}{2G}\end{bmatrix}$$

Solving for the stress tensor,


 * $$\left\{\sigma_{ij} \right\}_{6x1} = \begin{bmatrix}A^{-1}_{3x3} & 0_{3x3} \\ 0_{3x3} & B^{-1}_{3x3} \end{bmatrix}\left\{\varepsilon_{ij} \right\}_{6x1}$$

In order to verify this result, a new matrix is defined as


 * $$\begin{bmatrix}c\end{bmatrix}_{6x6} = \begin{bmatrix}A_{3x3} & 0_{3x3} \\ 0_{3x3} & B_{3x3} \end{bmatrix}$$ with inverse equal to $$\begin{bmatrix}c\end{bmatrix}^{-1}_{6x6} = \begin{bmatrix}A^{-1}_{3x3} & 0_{3x3} \\ 0_{3x3} & B^{-1}_{3x3} \end{bmatrix}$$

Where $$A^{-1}_{3x3} = \begin{bmatrix} \frac{E(\nu-1)}{(2\nu^2+\nu-1)} & \frac{-E\nu}{(2\nu^2+\nu-1)} & \frac{-E\nu}{(2\nu^2+\nu-1)}  \\ \frac{-E\nu}{(2\nu^2+\nu-1)} & \frac{E(\nu-1)}{(2\nu^2+\nu-1)}  & \frac{-E\nu}{(2\nu^2+\nu-1)}  \\ \frac{-E\nu}{(2\nu^2+\nu-1)} & \frac{-E\nu}{(2\nu^2+\nu-1)} & \frac{E(\nu-1)}{(2\nu^2+\nu-1)}\end{bmatrix}$$ and $$B^{-1}_{3x3} = \begin{bmatrix}2G & 0 & 0 \\ 0 & 2G & 0 \\ 0 & 0 & 2G\end{bmatrix}$$

Due to the definition of the inverse, it is known that the multiplication of matrix c with its inverse will yield the identity matrix, I.


 * $$\begin{bmatrix}c\end{bmatrix}_{6x6}\begin{bmatrix}c\end{bmatrix}^{-1}_{6x6} = \begin{bmatrix}I\end{bmatrix}_{6x6}$$


 * $$\begin{bmatrix}c\end{bmatrix}\begin{bmatrix}c\end{bmatrix}^{-1} = \begin{bmatrix}AA^{-1} & 0 \\ 0 & BB^{-1} \end{bmatrix} = \begin{bmatrix}I\end{bmatrix}$$

Now, the stress-strain relation can be found using the inverted matrices, A and B.


 * $$\left\{\sigma_{ij} \right\} = \begin{bmatrix}A^{-1} & 0 \\ 0 & B^{-1} \end{bmatrix}\left\{\varepsilon_{ij} \right\} = \begin{bmatrix} \frac{E(\nu-1)}{(2\nu^2+\nu-1)} & \frac{-E\nu}{(2\nu^2+\nu-1)} & \frac{-E\nu}{(2\nu^2+\nu-1)} & 0 & 0 & 0\\ \frac{-E\nu}{(2\nu^2+\nu-1)} & \frac{E(\nu-1)}{(2\nu^2+\nu-1)}  & \frac{-E\nu}{(2\nu^2+\nu-1)} & 0 & 0 & 0\\ \frac{-E\nu}{(2\nu^2+\nu-1)} & \frac{-E\nu}{(2\nu^2+\nu-1)} & \frac{E(\nu-1)}{(2\nu^2+\nu-1)}  & 0 & 0 & 0\\ 0 & 0 & 0 & 2G & 0 & 0\\ 0 & 0 & 0 & 0 & 2G & 0\\ 0 & 0 & 0 & 0 & 0 & 2G \end{bmatrix}\begin{Bmatrix}0 \\ 0 \\ 0 \\ 0 \\ \varepsilon_{31} \\ \varepsilon_{12} \end{Bmatrix} \;$$

Equation of Equilibrium Stresses
We now want to look at a non uniform stress field. To accomplish this we must first start with the basics and consider a 1-D case as a model. We will have to consider both a Uniform and Nonuniform case. The uniform case can be seen in Figure 1, while the nonuniform case can be seen in Figure 2. The uniform stress includes a single force acting on the end of a cantilever beam, while the nonuniform case has a varying shear stress acting horizontally across the top face of the beam along with the force on the end of the beam.



In the nonuniform case, the shear stress f(x) changes along the beam, thus f(x) is not a constant.

Bidirectional Bending


The recipe for Bidirectional Bending can be seen in Figure 3. The figure shows nonuniform stresses acting on a nonuniform cylinder. The equations below show the Moments about the y and z axis from $$ \displaystyle \sigma_{xx}$$. The Area over which the following integrals will be integrated are $$ dA = dy dz\;$$.


 * $$ M_y = \int_A z \displaystyle \sigma_{xx} \,dA $$


 * $$ M_z = \int_A y \displaystyle \sigma_{xx} \,dA $$

The next step is to look at the Moment of Inertia Tensors. First it is helpful to look at how the tensors can be transformed into both uniform and indicial notations.


 * $$ I_y = I_{yy} = I_{22}\;$$


 * $$ I_z = I_{zz} = I_{33}\;$$


 * $$ I_{yz} = I_{23}\;$$

Now from Figure 3, the Moment of Inertia Tensors can be determined.


 * $$ I_y = \int_A z^2 \,dA$$


 * $$ I_z = \int_A y^2 \,dA$$


 * $$ I_{yz} = \int_A yz \,dA$$

Now we will examine Hooke's Law:


 * $$\displaystyle \sigma_{xx} = E \epsilon_{xx}$$

Another way of looking of looking at Hooke's Law:


 * $$\displaystyle \sigma_{xx} = M_yy + M_zz$$


 * $$\displaystyle \sigma_{xx} = (\frac{I_y M_z - I_{yz} M_z}{I_y I_z - I_{yz}^2}) y + (\frac{I_z M_z - I_{yz} M_{z}}{I_y I_z - I_{yz}^2}) z$$

Note that the denominator is the determinant of the moment of inertia matrix.

Note that the neutral axis where $$\displaystyle \sigma_{xx} = 0$$. Thus Hooke's Law can be transformed into the following:


 * $$\displaystyle \sigma_{xx} = M_yy + M_zz = 0 $$


 * $$ z = \frac{-M_y}{M_z} y$$

By looking at the simple geometry of Figure 3:


 * $$ tan \displaystyle \beta = \frac{M_y}{M_z}$$

Thus the equation can be simplified:


 * $$ z = (tan\displaystyle \beta)y$$


 * {| class="collapsible collapsed"

!HW Problem - Equations (2.22) and (2.23) in Indicial Form
 * Equation 2.22 can be expressed as
 * Equation 2.22 can be expressed as


 * $$\frac{\partial\tau_{xy}}{\partial x} + \frac{\partial\sigma_{yy}}{\partial y} + \frac{\partial\tau_{zy}}{\partial z} = 0$$

and can be written in indicial form by replacing the subscripts


 * $$\;x = 1$$
 * $$\;y = 2$$
 * $$\;z = 3$$

And noting that $$\tau_{ij}$$ is really merely $$\sigma_{ij}$$ with a subscripts where $$i \neq j$$ since the first subscript is normal to the facet and the second subscript is the direction of the stress. If the subscript which is normal to the facet is not equal to the subscript which signifies the direction of the stress, then the force must be a shear force. Therefore, the notation using $$\tau$$ can be replaced by only using $$sigma$$. Using indicial form, the following is shown


 * $$\frac{\partial\sigma_{12}}{\partial x_1} + \frac{\partial\sigma_{22}}{\partial x_2} + \frac{\partial\sigma_{32}}{\partial x_3} = 0$$

Similarly, for equation 2.23


 * $$\frac{\partial\tau_{xz}}{\partial x} + \frac{\partial\tau_{yz}}{\partial y} + \frac{\partial\sigma_{zz}}{\partial z} = 0$$

can be reduced to


 * $$\frac{\partial\sigma_{13}}{\partial x_1} + \frac{\partial\sigma_{23}}{\partial x_2} + \frac{\partial\sigma_{33}}{\partial x_3} = 0$$


 * }


 * {| class="collapsible collapsed"

!HW Problem - General Form of Equations (2.22) and (2.23) in Indicial Form The indicial forms of equations 2.22 and 2.23 can be expressed in a general form using


 * $$\sum_{i=1}^3 \frac{\partial\sigma_{ij}}{\partial x_i} = 0$$

for
 * $$j =\; 1, 2, 3$$

For equations 2.22 and 2.23, $$j = 2$$ and $$j = 3$$, respectively. When summing them up, the following is obtained


 * $$\sum_{i=1}^3 \frac{\partial\sigma_{i2}}{\partial x_i} = \frac{\partial\sigma_{12}}{\partial x_1} + \frac{\partial\sigma_{22}}{\partial x_2} + \frac{\partial\sigma_{32}}{\partial x_3} = 0$$

and


 * $$\sum_{i=1}^3 \frac{\partial\sigma_{i3}}{\partial x_i} = \frac{\partial\sigma_{13}}{\partial x_1} + \frac{\partial\sigma_{23}}{\partial x_2} + \frac{\partial\sigma_{33}}{\partial x_3} = 0$$


 * }

Equation of Equilibrium
The goal of this section is to prove that the equation below must be equal to zero in order for the body to be in equilibrium. This equation is thus called the equilibrium equation. The equilibrium equation in the x-direction is:


 * $$\frac{\partial \sigma _{xx}}{\partial x}+ \frac{\partial \sigma _{yx}}{\partial y}+\frac{\partial \sigma _{zx}}{\partial z}=0$$

This equation can be expressed in indicial notation with $$x=1\;$$, $$y=2\;$$, and $$z=3\;$$ and compressed as follows:


 * $$\sum_{i=1}^{3}{\frac{\partial \sigma _{ij}}{\partial x_i}}=0$$ for $$j=1, 2, 3\;$$

We can also note that in the case of pure torsion, the normal stress component $${\sigma _{xx}}\;$$ is equal to zero, so the equilibrium equation simplifies to:


 * $$\frac{\partial \sigma _{yx}}{\partial y}+\frac{\partial \sigma _{zx}}{\partial z}=0$$

Before moving into the derivation of the equilibrium equation for the 3-D case, The 1-D case shown in Figure 4 should first be examined. In order for this stress element to be in equilibrium, the sum of the forces must equal zero.


 * $$\sum{F_x}=0=-\sigma (x)A+\sigma (x+dx)A+f(x)dx$$
 * $$0= A \left[ \sigma (x+dx)-\sigma (x)\right]+f(x)dx$$

Now, recall the Taylor Series expansion:


 * $$f(x+dx)=f(x) + \frac{df(x)}{dx}(dx) + \frac{1}{2}\frac{d^2f(x)}{dx^2}(dx)^2 + ...$$



We can replace the term in square brackets with the first order Taylor Series expansion to simplify the equation as follows:


 * $$\frac{d \sigma }{dx}+\frac{f(x)}{A}=0$$

In this equation, $$f(x)\;$$ is in units of force per length, and the term $$\frac{f(x)}{A}$$ is called the applied load. We will now extend this derivation into the three-dimensional case. We will consider a stress element in a nonuniform stress field, without an applied load, and focus only on the stress components in the x-direction to avoid over-complicating the problem. The derivation will yield the same results if we consider the stresses acting in the y- and z-directions. The 3-D problem is depicted in Figure 5. Note that the first subscript in the stress terms defines the direction normal to the plane of the stress, and the second subscript defines the direction of the stress. The origin of the coordinate system is at the location $$(x,y,x)\;$$.

The sum of the forces is expressed as:


 * $$\sum{F_x}=0=dydz\left[-\sigma _{xx}(x,y,z) + \sigma _{xx}(x+dx,y,z)\right] + dxdz\left[-\sigma _{yx}(x,y,z) + \sigma _{yx}(x,y+dy,z) \right] + dxdy\left[-\sigma _{zx}(x,y,z) + \sigma _{zx}(x,y,z+dz) \right]$$

Simplifying the equation using the Taylor Series expansion, as in the 1-D case, yields the following:


 * $$0=(dxdydz)\left[\frac{\partial \sigma _{xx}}{\partial x} + \frac{\partial \sigma _{yx}}{\partial y} + \frac{\partial \sigma _{zx}}{\partial z}\right]$$

The volume term cancels out, and we are left with the equilibrium equation for stresses in the x-direction for the 3-D problem:


 * $$0=\frac{\partial \sigma _{xx}}{\partial x} + \frac{\partial \sigma _{yx}}{\partial y} + \frac{\partial \sigma _{zx}}{\partial z}$$

MATLAB Airfoil Problem Continued
The MATLAB airfoil problem, continued for HW 5, consisted of two parts. Both parts dealt with analysis of the effects of bidirectional bending on an arbitrary 4-digit NACA airfoil. Part I refers to an idealized case where only the axial effects of the stringers are considered. Part 2 considers the "full-blown" case - the bending effects of the airfoil skin, spar webs at 1/4 and 3/4 chord, and the stringers are considered. Below is a list of given and researched values used during the problem:


 * $$ A_B = 2 x 10^{-4}\;m^{2}\;$$ - Cross sectional area of the stringer at point B
 * $$ A_E = 2 x 10^{-4}\;m^{2}\;$$ - Cross sectional area of the stringer at point E
 * $$ A_F = 1 x 10^{-4}\;m^{2}\;$$ - Cross sectional area of the stringer at point F
 * $$ A_H = 1 x 10^{-4}\;m^{2}\;$$ - Cross sectional area of the stringer at point H
 * $$ t_s = .002\;m\;$$ - Airfoil skin thickness
 * $$ t_w = .003\;m\;$$ - Spar web thickness
 * $$ t_{str} = .005\;m\;$$ - Stringer flange thickness
 * $$ M_y = -1250\;N\cdot m\;$$ - y-direction bending moment (compressive)
 * $$ M_z = 500\;N\cdot m\;$$ - z-direction bending moment (tensile)
 * $$ \sigma_{u} = 1860\;N/m\;$$ - Ultimate tensile stress for steel 300M

Part I: Ideal Case


Figure 6 shows a generic airfoil with spars located at 1/4 and 3/4 chord. The points B, E, H, and F, representing the intersection of the spars with the airfoil contour, are labeled. For the purposes of Part I, these 4 points represent the location of the idealized stringers. The stringer cross-section is not considered, and each stringer is represented as a "point area" during the calculations. Below is the full MATLAB code for this part followed by a sample output run.


 * {| class="collapsible collapsed"

!MATLAB Code Part I - Bending Effects, Only Considering Stringers
 * 
 * 
 * }


 * {| class="collapsible collapsed"

!MATLAB Code Part I Output >> vqcrew_4str(.5,2,4,15,1000);
 * 
 * 
 * }



Figuree 7 is the output plot from the MATLAB code for Part I. Several features are notable. First, the solid black circles located at the spar-airfoil intersections represent each stringer. The blue dashed lines represent the airfoil skin and spar webs. The red dash-dotted line running across the airfoil represents the neutral axis, and the black crosshair is the centroid of the stringers. The solid red circle represents the location of maximum normal stress on the airfoil.

A notable result from the output are the ultimate bending moments for the supplied ultimate stress. The calculated moments are of opposite sign to the original given moments. This is clearly explained by the result of the maximum bending stress analysis for the original loading condition. The determined stress was negative, indicating a compressive stress. The provided ultimate stress was given as a tensile stress, and since the neutral axis and thus the point of maximum normal stress do not change for this portion of the problem, the signs on the moments must switch to produce a tensile stress at that point. This also explains why much smaller magnitude moments produce a stress in excess of the ultimate stress by an order of magnitude - the ultimate stress provided is for tensile stresses, not compressive.

Part II: Full-Blown Case


For part II, the entirety of the airfoil was considered, including the cross-section of the stringers. In order to generate this cross-section, the assumption was made that the stringers were manufactured by stamping from a metal sheet of consistent thickness $$t_{str}\;$$. The following equation was used to determine the length of one flange of each stringer:


 * $$L_{brack} = .5*(A_{str}/2)/t_{str}\;$$

The stringers were assumed to be flush with the inner airfoil contour and the spar web wall. One stringer would consist of two "L-shaped" brackets on either side of the spar. Since the airfoil is assumed to have curvature, the brackets will not be identical, and the assumed manufacturing process means the angle between the two flanges will not be 90 degrees.

Figures 8 and 9 shows a generic depiction at both the 1/4 chord spar and the 3/4 chord spar. The diagram of the 1/4 chord spar is more representative of the general stringer shape. The stringers for the 3/4 chord spar are presented differently, however. The result of the above equation for the two stringers at the 3/4 spar was a bracket length of .005 m, which is identical to the provided stringer thickness. If one goes by the geometry provided in Figure 8, this essentially provides a quadrilateral cross-section, which is what is drawn in Figure 9 and also used in the code for Part II, which is provided below. After the code is a sample output run.


 * {| class="collapsible collapsed"

!MATLAB Code Part II - Full-Blown NACA Airfoil
 * 
 * 
 * }


 * {| class="collapsible collapsed"

!MATLAB Code Part II Output
 * 
 * 
 * }



Figure 10 depicts the full output plot of the Part II code. The airfoil skin and spar webs were plotted as solid blue lines. The stringer cross-sections are plotted as small black circles, with each portion of the stringer wall represented with 5 points for ease of viewing. Zoomed in plots of the 1/4 spar and 3/4 spar are provided in Figures 11 and 12 respectively. The neutral axis is again plotted as a red dash-dotted line.

There are two centroids plotted on this airfoil; a zoomed in view of them is available in Figure 13. The green crosshair represents the centroid of the stringers when the stringer cross-sections are considered. The blue crosshair is the crosshair calculated in Part I. Also plotted is a red solid circle again representing the location of maximum bending stress on the airfoil, and a green solid circle represents the location of maximum bending stress in a stringer. Comparison and conclusions about these results follows in the next section.

Comparison and Conclusions
Below is a table showing selected data from the output of both cases. The stringer centroids, locations of maximum bending stresses, and neutral axis angle for both cases are more or less identical, which is not surprising given both are based off of the same airfoil with the same input data. The most apparent differences are in the data for ultimate bending moments and maximum normal bending stress. The higher maximum stress at lower supplied moment indicates the full-blown airfoil is slightly less efficient at carrying bending loads than the idealized airfoil. This is probably due primarily to the stringers being much, much better at carrying bending loads than either the spars or the airfoil skin. The only other notable point from the analysis is the stringer which carries the most stress in the full-blown case, which was found to be stringer B. This is not surprising, as this stringer is closest to the position of maximum thickness, which in combination with the positive slope of the neutral axis means it would be the farthest stringer from the neutral axis, and thus experience the highest stress.


 * {| border=1


 * bgcolor="red"|
 * bgcolor="#99ccff" style="text-align: center;" |Ideal Case
 * bgcolor="#99ccff" style="text-align: center;" |Full-Blown Case
 * bgcolor="red"|Stringer centroid (meters,meters)
 * bgcolor="#99ccff" style="text-align: center;" | (0.2083 0.0080)
 * bgcolor="#99ccff" style="text-align: center;" | (0.2088 0.0080)
 * bgcolor="red"|Moment of Inertia Tensor Components (m^4,m^4,m^4)
 * bgcolor="#99ccff" style="text-align: center;" | (6.6875e-007,3.4375e-005,9.3322e-007)
 * bgcolor="#99ccff" style="text-align: center;" | (6.2884e-007,1.4533e-004,4.8673e-006)
 * bgcolor="red"|Neutral Axis Slope (radians)
 * bgcolor="#99ccff" style="text-align: center;" | 0 .0345
 * bgcolor="#99ccff" style="text-align: center;" | 0.0347
 * bgcolor="red"|Maximum normal bending stress (Pascals)
 * bgcolor="#99ccff" style="text-align: center;" | -8.1995e+007 Pascals
 * bgcolor="#99ccff" style="text-align: center;" | -1.1641e+008 Pascals
 * bgcolor="red"|Location of maximum normal bending stress (meters,meters)
 * bgcolor="#99ccff" style="text-align: center;" | (0.1436 0.0467)
 * bgcolor="#99ccff" style="text-align: center;" | (0.1434 0.0478)
 * bgcolor="red"|Ultimate Bending Moments (N-m, N-m)
 * bgcolor="#99ccff" style="text-align: center;" | (28355,-11342)
 * bgcolor="#99ccff" style="text-align: center;" | (20970,-8388)
 * }
 * bgcolor="#99ccff" style="text-align: center;" | (0.1434 0.0478)
 * bgcolor="red"|Ultimate Bending Moments (N-m, N-m)
 * bgcolor="#99ccff" style="text-align: center;" | (28355,-11342)
 * bgcolor="#99ccff" style="text-align: center;" | (20970,-8388)
 * }
 * }
 * }

Contributing Team Members
Steven Hepsworth Eas4200c.f08.vqcrew.f 13:44, 6 November 2008 (UTC) Darin Toscano Eas4200c.f08.vqcrew.d 22:52, 6 November 2008 (UTC) Kevin Klauk Eas4200c.f08.vqcrew.E 16:55, 7 November 2008 (UTC) Adam Edstrand Eas4200c.f08.vqcrew.a 18:12, 7 November 2008 (UTC) John Saxon Eas4200c.f08.vqcrew.c 19:59, 7 November 2008 (UTC) Keith Javier Stober Eas4200c.f08.vqcrew.b 16:38, 8 December 2008 (UTC)