User:Eas4200c.f08.aeris.krammer/HW4

 Comment: This page is identical to the previous version, except that notes on the Strain Matrix have been added (Item 7 in the Contents tab). Here you may find a comparison between this and the original version. Eas4200c.f08.aeris.mcrae. 15:44, 29 October 2008 (UTC)

Multicell Airfoil
Let's take a specific example first to find $$\displaystyle \theta$$ as a function of $$\displaystyle T$$ and $$\displaystyle J$$ (torsional constant) First, $$T=T_{1}+T_{2}=2q_{1}\bar{A_{1}}+2q_{2}\bar{A_{2}}  (1)$$ Where, $$\bar{A_{1}}=ac$$ $$\bar{A_{2}}=bc$$ Then, $$\theta_{1}=\frac{1}{2G\bar{A_{1}}}\oint_{}^{}{}\frac{q_{1}}{t_{1}(s)}ds=\frac{1}{2G\bar{A_{1}}}\left[\frac{2q_{1}b}{t_{1}} +\frac{q_{1}c}{t_{1}}+ \frac{(q_{1}-q_{2})c}{t_{1,2}}\right]   (2)$$ $$\theta_{2}=\frac{1}{2G\bar{A_{2}}}\oint_{}^{}{}\frac{q_{2}}{t_{2}(s)}ds=\frac{1}{2G\bar{A_{2}}}\left[\frac{2q_{2}b}{t_{2}} +\frac{q_{2}c}{t_{2}}+ \frac{(q_{2}-q_{1})c}{t_{1,2}}\right]   (3)$$ $$\displaystyle C_{1}$$ and $$\displaystyle C_{2}$$ have the same rate of twist angle, therefore $$\displaystyle \theta _{1}=\theta _{2}$$ Using equation (1), we can find an expression for $$\displaystyle q_{1}$$ and $$\displaystyle q_{2}$$ in terms of $$\displaystyle T$$: $$\displaystyle q_{1}=\beta _{1}T$$ $$\displaystyle q_{2}=\beta _{2}T$$ Where $$\displaystyle \beta _{1}$$ and $$\displaystyle \beta _{2}$$ is some number. Next, we can use equations (2) and (3) to find the to find the expression for $$\displaystyle \theta$$ and $$\displaystyle T$$, $$\displaystyle \theta =\theta _{1}=\theta _{2}=\frac{T}{2GJ}$$, and deduce $$\displaystyle J$$. Calculating with numerical values for the above equations: $$\displaystyle A_{1}=0.3m * 0.4 m =0.12m^2$$ $$\displaystyle A_{2}=0.6m * 0.4 cm =0.24m^2$$ Therefore,  $$\displaystyle T=0.24q_{1}+0.48q_{2}$$ Assuming $$\displaystyle G_{1}=G_{2}=G$$, we can solve for $$\displaystyle \theta _{1}$$ and $$\displaystyle \theta _{2}$$ $$ \theta _{1}=\frac{1}{0.24G}[\frac{1.2q_{1}}{0.003}+\frac{0.4q_{1}}{0.003}+\frac{0.4(q_{1}-q_{2})}{0.004}]$$ $$ \theta _{2}=\frac{1}{0.48G}[\frac{1.2q_{2}}{0.005}+\frac{0.4q_{2}}{0.005}+\frac{0.4(q_{2}-q_{1})}{0.004}]$$ Since $$\displaystyle \theta _{1}=\theta _{2}$$ equating the two equations above yields to $$\displaystyle 1366.6q_{1}=620q_{2}$$ Using equation $$\displaystyle (1)$$, we can express $$\displaystyle q_{1}$$ and $$\displaystyle q_{2}$$ as a function of $$\displaystyle T$$ $$\displaystyle q_{1}=0.77T$$  $$ \displaystyle q_{2}=1.698T$$ Where $$\displaystyle 0.77=\beta _{1}$$ and $$\displaystyle 1.698=\beta _{2}$$

Derivation of the Twist Angle in the Case of a Uniform, Non-Circular Cross Sectional Bar with Varying Shear Flow
Regarding the previous lecture, once $$\theta$$ is found in terms of T, use $$T=GJ\theta$$ or $$J=\frac{T}{G\theta}$$ to find the torsional constant J.  For this homework assignment, the NACA 2415 airfoil, referring to lecture 16.2, is partitioned into three cells. The cells are separated by two walls of assigned thickness. To perform the Matlab analysis similar to the homework 3, more information must be provided in order to properly calculate the "average" areas $$\bar{A}_1, \bar{A}_2, and \bar{A}_3$$ corresponding to the three cells of the airfoil and along with the torque (not indicated but will be treated as a constant), thickness of the walls and shell, contour length and shear modulus will allow us to calculate the shear flow on the airfoil. First, one wall is .25 times the chord length away from the leading edge, and a second wall is .75 times the chord length away from the leading edge. The shear modulus is not needed because the equations used will cancel out these terms. There will be three unknowns $$q_1,q_2, and q_3$$ and using the three following equations, these unknowns will be solved for. $$T=2\sum_{i=1}^3{q_i\bar{A}_i}$$ $$\displaystyle \theta_1=\theta_2$$ $$ \displaystyle \theta_2=\theta_3$$ Then, the torsional constant J will be found by the equation mentioned above. As we have already finished deriving the follow equations, $$T=GJ\theta, T=2q\bar{A}$$ Now we will derive the twist angle equation for a uniform, non-circular cross sectional with an engineering (ad hoc) method. $$\theta=\frac{1}{2G\bar{A}}\oint\frac{q}{t}ds$$ Let's consider this uniform, non-circular cross sectional bar subject to a torque. The bottom two pictures describe this type of bar with infinitesimally small axial length dx and the cross section is shown in the (y,z) plane. <p style="text-align:center;"> Displacement PP' is due to alpha, and approximating alpha as very small, we have the following expression for PP' and OP in a simple geometric triangular relation. <p style="text-align:center;">$$\frac{PP'}{OP}=tan(\alpha)$$ Now, project the displacement PP' in the direction orthogonal to OP'. In other words, <p style="text-align:center;">$$\displaystyle PP' = PP'' cos( \alpha)$$ $$\displaystyle PP=(OPtan(\alpha))cos(\alpha)=(OPcos(\alpha))tan(\alpha)=OPtan(\alpha)$$ Recall, that OP=R, the radial coordinate of the point P and thus $$\rho$$ must be the perpendicular distance from the origin to P which equals, in this case, OP. Knowing that, the following expression can be written. <p style="text-align:center;">$$\displaystyle PP''=(Rcos(\alpha))tan(\alpha)=\rho \alpha$$ Here, PP" is the displacement of P in the direction "tangent" to the lateral surface of the bar. Strain is $$\gamma=\frac{PP}{dx}$$ where dx is the distance between very close planes on the bar. <p style="text-align:center;">$$\gamma=\frac{PP}{dx}=\frac{\rho \alpha}{dx}=\rho \theta$$ Here theta is the rate of twist angle.

Formal Justification by Elasticity Theory
<p style="text-align:center;">$$\displaystyle \tau =G\gamma=G\rho \theta$$ Where the rate of twist $$\displaystyle \theta =\frac{d\alpha }{dx}$$ Integrating along the contour wall <p style="text-align:center;">$$\oint_{C}^{}{}\tau (s) ds = G\theta (x)\oint_{C}^{}{}\varrho (s) ds$$ Where $$\displaystyle \tau (s) ds$$ is $$\displaystyle \frac{q(s)}{t(s)}$$ and $$\oint_{C}^{}{}\varrho (s) ds$$ is  $$\displaystyle 2\bar{A}$$ Hence we get the following expression for $$\displaystyle \theta$$ <p style="text-align:center;">$$\theta =\frac{1}{2GA}\oint_{}^{}{}\frac{q}{t}ds$$ What is ad hoc about this derivation and the derivation of $$\displaystyle T=2q\bar{A}$$? 1) Strain must be obtained using the displacement of $$\displaystyle \vec{P}$$ in the direction tangent to $$\displaystyle C$$ at $$\displaystyle \vec{P}$$ but $$\displaystyle \vec{PP'}$$ is not necessarily tangent to $$\displaystyle C$$, although it is really close. 2) $$\displaystyle \tau =\frac{q}{s}$$ obtained from the ad-hoc assumption that $$\displaystyle \tau$$ was uniform across the wall thickness. Thus, we have the formal derivation by elasticity theory.

Roadmap: Kinematic Assumption
The kinematic assumptions are: <p style="text-align:center;">$$\displaystyle u_{x}(y,z)=\theta \psi (y,z)$$ (Where $$\displaystyle \theta$$ is considered to be constant with respect to $$\displaystyle x$$) <p style="text-align:center;">$$\displaystyle u_{y}(x,z)=-\theta xz$$ <p style="text-align:center;">$$\displaystyle u_{z}(x,y)=\theta xy$$ The following relationship can be established: <p style="text-align:center;">$$\displaystyle \epsilon _{xx}=\epsilon _{yy}=\epsilon _{zz}=\gamma _{yz}=0$$ This relation can be proved by taking partial derivatives: <p style="text-align:center;">$$\epsilon _{xx}=\frac{\partial u_{x}(y,z)}{\partial x}=0$$ <p style="text-align:center;">$$\epsilon _{yy}=\frac{\partial u_{y}(x,z)}{\partial y}=0$$ <p style="text-align:center;">$$\epsilon _{zz}=\frac{\partial u_{z}(x,y)}{\partial z}=0$$ <p style="text-align:center;">$$\gamma _{yz}=\frac{\partial u_{y}(x,z)}{\partial z}+\frac{\partial u_{z}(x,y)}{\partial y}=0$$ (Where $$\frac{\partial u_{y}(x,z)}{\partial z}=-\theta x$$ and $$\frac{\partial u_{z}(x,y)}{\partial y}=+\theta x$$)

NACA Airfoil: Triangle Rule
The triangle rule can be used to find the area of any 4-digit NACA Airfoil from an observer at any point. The total area would be found by:

$$A = A_{1} + A_{2}$$ $$ = A_{1} - |A_{2}|$$

To find the areas of $$A_{1}$$ and $$A_{2}$$, the triangle rule is used by the following:

This method works effectively because the observer is noted as the origin and each triangle stems from that origin. This allows for infinitely many triangles to be measured from a single point and an area to be calculated correctly and accurately. Another advantage to the triangle rule is that it can be used in any case.

 EXAMPLE 1 

Although it may be easier to integrate a rigid body using other methods (i.e. Trapezoidal Rule), curvature changes, like the one shown above, may not be very easy to compute. The curve causes for tricky calculations which can be avoided by using triangles.

Single Cell Airfoil


Problem Statement: If shear flow is assumed to be constant (i.e. $$q = q_{1} = q_{2} = q_{3}$$), what is the angle of twist and the allowable torsion, $$T_{allow}$$.

First, we must calculate the rate of twist, θ,  θ = (q/2GA)Σ lj/tj . The variables are: 1. q = shear flow 2. lj = length 3. G = shear modulus of cell 4. tj = thickness in cell (curves along cell wall) 5. A = cell area

Applied to this single cell, the equation is:  θ = (q/2GA)[π*(b/2)/t1 + a/t2 + sqrt(a2 + b2)/t3] 

Once the rate of twist is calculated we can then find the maximum shear stress, τmax. The formula for calculating τmax is:  τmax = q / min{t1, t2, t3}  where min{t1, t2, t3} refers to the smallest thickness value of the the three used.

Because this is a single cell,  τmax = τallow  and since $$q = T/2A$$, then  Tallow = 2A*τallow*min{t1, t2, t3} 

Strain Matrix
In a 3-D body under strain, this strain can be broken down into a component matrix, as seen below.

$$\varepsilon =\begin{bmatrix} \varepsilon _{xx} & \varepsilon _{xy} & \varepsilon _{xz}\\ \varepsilon _{yx} & \varepsilon _{yy} & \varepsilon _{yz}\\ \varepsilon _{zx} & \varepsilon _{zy} & \varepsilon _{zz} \end{bmatrix} =\begin{bmatrix} \varepsilon _{11} & \varepsilon _{12} & \varepsilon _{13}\\ \varepsilon _{21} & \varepsilon _{22} & \varepsilon _{23}\\ \varepsilon _{31} & \varepsilon _{32} & \varepsilon _{33} \end{bmatrix} =\left[\varepsilon _{ij} \right]$$

Due to symmetry, there are a several equalities we can make.

$$\varepsilon _{xy}=\varepsilon _{yx}$$

$$\varepsilon _{yz}=\varepsilon _{zy}$$

$$\varepsilon _{xz}=\varepsilon _{zx}$$

We can now see that this 3x3 strain matrix has only 6 components, $$\varepsilon _{xx}, \varepsilon _{yy},, \varepsilon _{zz}, \varepsilon _{xy}, \varepsilon _{yz}, \varepsilon _{xz}$$.

The equation for each component can be expressed as $$\varepsilon _{ij}=\frac{1}{2}\left(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\right)$$

Stress Tensor Matrix

Similarly, we can make some equalities to the stress tensor matrix, which is seen below.

$$\sigma =\begin{bmatrix} \sigma _{xx} & \sigma _{xy} & \sigma _{xz}\\ \sigma _{yx} & \sigma _{yy} & \sigma _{yz}\\ \sigma _{zx} & \sigma _{zy} & \sigma _{zz} \end{bmatrix} =\begin{bmatrix} \sigma _{11} & \sigma _{12} & \sigma _{13}\\ \sigma _{21} & \sigma _{22} & \sigma _{23}\\ \sigma _{31} & \sigma _{32} & \sigma _{33} \end{bmatrix} =\left[\sigma _{ij} \right]$$

$$\sigma _{xy}=\sigma _{yx}$$

$$\sigma _{yz}=\sigma _{zy}$$

$$\sigma _{xz}=\sigma _{zx}$$

We see that the stress tensor matrix, like the strain matrix, also has 6 components.

Stress-Strain Relation
Normal Strains: <p style="text-align:center;"> $$\epsilon _{xx}=\frac{\sigma _{xx}}{E}-\frac{\nu \sigma _{yy}}{E}-\frac{\nu \sigma _{zz}}{E}$$ <p style="text-align:center;"> $$\epsilon _{yy}=\frac{\sigma _{yy}}{E}-\frac{\nu \sigma _{xx}}{E}-\frac{\nu \sigma _{zz}}{E}$$ <p style="text-align:center;"> $$\epsilon _{zz}=\frac{\sigma _{zz}}{E}-\frac{\nu \sigma _{xx}}{E}-\frac{\nu \sigma _{yy}}{E}$$ Shear Strains: <p style="text-align:center;"> $$\gamma _{xy}=2\epsilon _{xy}=\frac{\tau _{xy}}{G}$$ <p style="text-align:center;"> $$\gamma _{yz}=2\epsilon _{yz}=\frac{\tau _{yz}}{G}$$ <p style="text-align:center;"> $$\gamma _{zx}=2\epsilon _{zx}=\frac{\tau _{zx}}{G}$$ Another form, voigt notation, is used to describe normal strain and shear strain in a column matrix. This allows $$\epsilon=\begin{bmatrix} \epsilon_{ij} \end{bmatrix}_{3\times 3}$$ and $$\sigma=\begin{bmatrix} \sigma_{ij} \end{bmatrix}_{3\times 3}$$ to be written as: <p style="text-align:center;"> $$\begin{Bmatrix} \epsilon _{ij} \end{Bmatrix}_{6\times 1}=\begin{Bmatrix} \epsilon _{11}\\ \epsilon _{22}\\ \epsilon _{33}\\ \epsilon _{23}\\ \epsilon _{31}\\ \epsilon _{12}\\ \end{Bmatrix}_{6\times 1} \begin{Bmatrix} \sigma _{ij} \end{Bmatrix}_{6\times 1}=\begin{Bmatrix} \sigma _{11}\\ \sigma _{22}\\ \sigma _{33}\\ \sigma _{23}\\ \sigma _{31}\\ \sigma _{12}\\ \end{Bmatrix}_{6\times 1}$$

Hooke's Law for Isotropic Elasticity: <p style="text-align:center;"> $$\begin{Bmatrix} \epsilon _{11}\\ \epsilon _{22}\\ \epsilon _{33}\\ \epsilon _{23}\\ \epsilon _{31}\\ \epsilon _{12}\\ \end{Bmatrix}_{6\times 1}=\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}_{6\times 6}\begin{Bmatrix} \sigma _{11}\\ \sigma _{22}\\ \sigma _{33}\\ \sigma _{23}\\ \sigma _{31}\\ \sigma _{12}\\ \end{Bmatrix}_{6\times 1} $$ which can also be written in terms of gamma and tau: <p style="text-align:center;"> $$\begin{Bmatrix} \epsilon _{11}\\ \epsilon _{22}\\ \epsilon _{33}\\ \gamma _{23}\\ \gamma _{31}\\ \gamma _{12}\\ \end{Bmatrix}_{6\times 1}=\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}_{6\times 6}\begin{Bmatrix} \sigma _{11}\\ \sigma _{22}\\ \sigma _{33}\\ \tau _{23}\\ \tau _{31}\\ \tau _{12}\\ \end{Bmatrix}_{6\times 1}$$

HW: Airfoil Matlab Code
For this homework assignment, our task was to modify the airfoil code from the last homework to make it work for multi-cell airfoils. The airfoil was to be divided at a quarter chord ($$c/4$$) and three-quarter chord ($$3c/4$$). The profile of the airfoil from last homework is shown below.

$$\bar{A}=0.0256 m^2$$ For this week we split the airfoil into three cells. The new profile for the airfoil is shown below.

$$\bar{A_{1}}=0.0072 m^2$$ $$\bar{A_{2}}=0.0159 m^2$$ $$\bar{A_{3}}=0.0026 m^2$$ The code had to be modified to calculate the area for each cell. To calculate the area for cell 1, we first found the closest location of $$y=c/4$$. Using the cross-product of two position vectors we were able to calculate the area to the leading edge of the airfoil. Next we calculated the area for cell 3 using a similar procedure. After finding the point of $$y=3c/4$$ we were able to calculate the area to the trailing edge. For the single-cell structure we obtain a torsional constant of 5.1118e-005. By adding two spars at c/4 and 3c/4 the torsional constant increases to 5.6071e-005. The additional spars did not increase the torsional constant by much.

Matlab Code Certification
I, the undersigned, certify that I can read, understand, and write matlab codes, and can thus contribute effectively to my team. Anett Krammer Eas4200c.f08.aeris.krammer 20:25, 23 October 2008 (UTC) Jin Yu Guan Eas4200c.f08.aeris.guan 20:29, 23 October 2008 (UTC) 02:45, 24 October 2008 (UTC) Nelson B Caceres Eas4200c.f08.aeris.caceres 14:13, 24 October 2008 (UTC) Josh Holladay Eas4200c.f08.aeris.holladay 19:36, 24 October 2008 (UTC) Jesse W. McRae Eas4200c.f08.aeris.mcrae 20:36, 28 October 2008 (UTC)

Contributions
Jin Yu Guan Eas4200c.f08.aeris.guan 19:01, 23 October 2008 (UTC) Anett Krammer Eas4200c.f08.aeris.krammer 19:55, 23 October 2008 (UTC) 02:44, 24 October 2008 (UTC) Nelson B Caceres Eas4200c.f08.aeris.caceres 14:13, 24 October 2008 (UTC) Josh Holladay Eas4200c.f08.aeris.holladay 19:36, 24 October 2008 (UTC) Jesse W. McRae Eas4200c.f08.aeris.mcrae 15:44, 29 October 2008 (UTC)