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

 Changes: The two MATLAB codes for calculation of the torsional constant J for both a single-celled and multi-celled airfoil have been modified and updated. The results of running the new codes have also been changed in the section.

Link to the pages of differences: Changes

Eas4200c.f08.vqcrew.c 15:49, 30 October 2008 (UTC)

Quadrature of a NACA Airfoil


The area of any shape can be determined by using a collection of many small triangles. The triangles are defined by an arbitrary origin, O, and an infinitesimal length on the shape. Since the infinitesimal length and perpendicular distance from the origin to the surface are known, the infinitesimal area, dA, can be calculated. Figure 1 shows this process for an airfoil.


 * $$ \bar{A}_{airfoil} = \bar{A}_{1} + \bar{A}_{2} = \bar{A}_{1} - \left|\bar{A}_{2}\right|$$

[[Image:vqcrew_triangle.png|thumb|right|200px|Fig. 2 - Determination of the infinitesimal area, dA, using the triangle method. ] In this case, the upper surface yields a positive area and the lower surface yields a negative area due to the location of the origin and the nature of the right-hand rule for vectors. The total area of the airfoil is computed as a sum of all of the small areas for each length along the airfoil contour.  Figure 2 provides a more detailed picture of how each area is calculated.  The example shown is for the lower surface, therefore the cross-product of the vectors  $$\textstyle \vec{r}$$  and $$\textstyle \vec{PQ}$$ is negative.  The negative cross-product leads to a negative area for the lower surface in this configuration.


 * $$\displaystyle d\vec{A} = \frac{1}{2}(\vec{r}\times\vec{PQ})$$



Problems with Trapezoidal Method of Quadrature
The process of dividing the area of a shape in to a series of small trapezoids may be an accurate way to determine the total area for some shapes in some configurations, but it is not a general solution. For example, Figure 3 shows the area of an airfoil being divided into trapezoids horizontally. For this shape in this orientation, the result will not be very accurate due to the errors near the trailing edge. The thin trailing edge is not conducive to being split up in this manner. It is true that if the trapezoids were aligned vertically the result would more accurately represent the true area, but a more elegant and beautiful solution can be reached. Using the triangle method described above, a more accurate result for the area can be found regardless of the location of the origin or orientation of the airfoil.

Single Cell Airfoil


For a single cell airfoil shown in Figure 4 the shear flow is constant.


 * $$\displaystyle q = q_{1} = q_{2} = q_{3}$$

Multi-cell Airfoil
In lecture 19, the professor introduced a specific problem in which he stated he would further generalize after the problem was solved. The image presented in Figure 5 shows a multi-cell airfoil which shows the shear flow through the airfoil as well as the various dimensions of the cells. The values given by the problem statement were as follows:




 * $$a \;= \;30$$ cm
 * $$b \;= \;60$$ cm
 * $$c \;= \;40$$ cm
 * $$t_1 \;= \;0.3$$ cm
 * $$t_2\; = \;0.5$$ cm
 * $$t_{12}\; = \;0.4$$ cm

The overall goal of the problem is to find $$\theta$$ as a function of $$T$$ and to deduce $$J$$ from that.


 * From lecture 17, it was shown that


 * $$T = 2\ \sum_{i=1}^{N_c} q_i\ \bar A$$

which resulted in what will further be defined as equation (1)


 * $$T = T_1 + T_2 = 2\ (q_1\ \bar A_1 + q_2\ \bar A_2)$$

where


 * $$\bar A_1 = a\ c$$
 * $$\bar A_2 = b\ c$$

It can be stated that the two cells have the same rate of twist angles, $$\theta\;$$. That being said, equation (2) can be shown that


 * $$\theta_1 = \theta_2 = ... = \theta_n\;$$

From lecture 17, page 2 of the notes, it can be shown that


 * $$\theta_i = \frac{1}{2\ G\ \bar A_i}\ \oint \frac{q_i}{t_i}\ ds$$

When substituting the values for the two cells, the following is obtained


 * $$\theta_1 = \frac{1}{2\ G\ \bar A_1}\ (\frac{2\ q_1\ a}{t_1} + \frac{q_1\ c}{t_1} + \frac{(q_1 - q_2)\ c}{t_{12}})$$
 * $$\theta_2 = \frac{1}{2\ G\ \bar A_2}\ (\frac{2\ q_2\ a}{t_2} + \frac{q_2\ c}{t_2} + \frac{(q_2 - q_1)\ c}{t_{12}})$$

Which will further be referred to as equation (3) and equation (4), respectively.

Now using equation (3) and equation (4), and using the relation


 * $$T = G\ J\ \theta$$

From this, $$J$$ can be deduced, resulting in


 * $$J = \frac{T}{G\ \theta}$$

Ad-hoc (Engineering) Derivation of Rate of Twist Formula
To derive the twist angle formula
 * $$\theta = \frac{1}{2G\bar{A}}\oint_{C}^{}{\frac{q}{t}ds}$$

we will consider a uniform bar with non-circular cross section subjected to twisting, as shown in Figure 6

To begin, the displacement $$PP'$$ is found to be


 * $$PP'=OPtan \mathbf{\alpha}$$

Then, project $$PP'$$ into the direction perpendicular to $$OP'$$ :
 * $$PP''=PP'cos \mathbf{\alpha}$$

Now the displacement of $$P$$ in the direction tangent to the lateral surface of the beam, $$PP$$, can be expressed in terms of $$OP$$ :
 * $$PP=(OPtan \mathbf{\alpha})cos \mathbf{\alpha} = OPcos \mathbf{\alpha}tan \mathbf{\alpha} = OPtan \mathbf{\alpha}$$

As in Figure 6, $$OP''$$ is defined as the perpendicular distance to the surface, $$\mathbf{\rho}$$. If we assume that the twist angle, $$\mathbf{\alpha}$$, is a small angle, $$PP''$$ simplifies to:
 * $$PP''=\mathbf{\rho}\mathbf{\alpha}$$

The shear strain, $$\mathbf{\gamma}$$, can be defined as the ratio of this displacement distance, $$PP''$$, to a small displacement in the axial (x) direction:
 * $$\gamma = \frac{PP''}{dx}$$

We can then substitute for $$PP''$$, and since we assumed $$\mathbf{\alpha}$$ is very small it can be replaced by it's differential element $$d\mathbf{\alpha}$$. We can also define the rate of twist, $$\mathbf{\theta}$$, as the differential twist angle $$d\mathbf{\alpha}$$ per differential change in the axial direction $$dx$$ :
 * $$\gamma = \frac{\rho d\alpha}{dx} = \rho \theta$$

We can use Hooke's Law to relate the shear stress to the shear strain:
 * $$\mathbf{\tau =G \gamma =G \rho \theta}$$
 * $$\mathbf{\tau(s) =G \rho (s) \theta (x)}$$

If we take the contour integral along C, the equation becomes:
 * $$\oint_{C}^{}{\tau(s)ds} = G \theta(x)\oint_{C}^{}{\rho (s)ds}$$

Next we can substitute in values from previous derivations to arrive at the final form for $$\mathbf{\theta}$$:
 * $$\oint_{C}^{}{\frac{q(s)}{t(s)}ds} = G \theta(x)2\bar{A}$$
 * $$\Rightarrow \theta = \frac{1}{2G\bar{A}}\oint_{C}^{}{\frac{q}{t}ds}$$

What makes this derivation "ad hoc"?

 * 1) The shear strain $$\mathbf{\gamma}$$ must be obtained using the displacement of $$P$$ in the direction tangent to $$C$$ at $$P$$, but $$PP''$$ is not necessarily tangent to C (but is close).
 * 2) The equation for the shear stress, $$\mathbf{\tau} =\frac{q}{t}$$, was obtained from an ad hoc assumption that $$\mathbf{\tau}$$ is uniform across the wall thickness.
 * 3) There is an inconsistency in the assumption on the size of $$\mathbf{\alpha}$$. To get line $$PP'$$, we assumed $$\mathbf{\alpha}$$ was small, but to get $$PP''$$ we assumed a finite $$\mathbf{\alpha}$$ to use the relation $$\mathbf{\rho}=OPcos \mathbf{\alpha}$$. We then reintroduced the small angle approximation for $$\mathbf{\alpha}$$ in the rest of the derivation.

Elastic Theory
For proper understanding of the elastic theory it can help to go back to the Kinematic Assumptions within the road map laid out in previous lectures and stated in the HW 3 report. Recall that in these assumptions that $$\displaystyle \theta$$ is considered to remain constant with respect to 'x' for a uniform bar.

A) Kinematic Assumptions


 * (1) $$u_x (y,z) = \displaystyle \theta \displaystyle \psi (y,z)$$


 * (2) $$u_y (x,z) = -\displaystyle \theta x z$$


 * (3) $$u_z (x,y) = +\displaystyle \theta x y$$



These equations vary slightly from the equations used in the textbook by C.T. Sun, this is due to the use of a different coordinated system. To transform the textbook equations to the uniform notation equations used above it is advised to use cyclic permutation, as seen in Figure 7. The basic relationship can be followed by the figure, but the arrows are used to represent the conversion between the x-axis of C.T. Sun becomes the y-axis in uniform notation and likewise for the other axis.

Other important relations can be generated from use of the Kinematic Assumptions and that is the relationship to certain strains and strain rate, giving the following results:


 * $$\displaystyle \epsilon_{xx} = \displaystyle \epsilon_{yy} = \displaystyle \epsilon_{zz} = \displaystyle \gamma_{yz} = 0$$

Formal Derivation
Now looking at the stains in 3-D. There are a total of six strain components for a 3 dimensional object due to symmetry.


 * $$\mathbf{\displaystyle \epsilon} = \begin{bmatrix}

\epsilon_{xx} & \epsilon_{xy} & \epsilon_{xz} \\ \epsilon_{yx} & \epsilon_{yy} & \epsilon_{yz} \\ \epsilon_{zx} & \epsilon_{zy} & \epsilon_{zz} \end{bmatrix}$$

Now that we have a 3x3 matrix with 9 coefficients, we need to change the notation to that of indicial notation.


 * $$ x \Leftrightarrow 1, y \Leftrightarrow 2 , z \Leftrightarrow 3 $$

Now we can write the matrix in indicial notation.


 * $$\mathbf{\displaystyle \epsilon} = \begin{bmatrix}

\epsilon_{11} & \epsilon_{12} & \epsilon_{13} \\ \epsilon_{21} & \epsilon_{22} & \epsilon_{23} \\ \epsilon_{31} & \epsilon_{32} & \epsilon_{33} \end{bmatrix}$$ $$ = \displaystyle \epsilon_{ij}$$

$$\displaystyle \epsilon_{ij}$$ is known at tensorial notation. For this matrix i,j range from 1 to 3, where i is the row number in the matrix and j in the column designation.

Now that the matrix for $$\displaystyle \epsilon_{ij}$$ is formed, the symmetry of the matrix can now be calculated.


 * $$\displaystyle \epsilon_{ij} = \frac{1}{2} (\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i})$$

The coordinates x,y,z can now be rewritten: $$ x \Leftrightarrow x_1, y \Leftrightarrow x_2 , z \Leftrightarrow x_3 $$


 * $$\epsilon_{11} = \displaystyle \epsilon_{xx} = \frac{1}{2} (\frac{\partial u_1}{\partial x_1} + \frac{\partial u_1}{\partial x_1})$$


 * $$\displaystyle \epsilon_{xx} = \frac{\partial u_1}{\partial x_1} \Rightarrow \frac{\partial u_x}{\partial x}$$



Stress Tensor
The stress tensor can be looked at much like the the strain relations as a 3x3 matrix. Also by symmetry there will be only 6 independent components of $$\displaystyle \sigma$$ in 3-D.


 * $$\displaystyle \sigma = [\sigma]_{3x3}$$

Based of the relations between the stress and strain, several stresses and a shear strain can be found by kinematic assumptions.


 * $$\displaystyle \sigma_{xx} = \displaystyle \sigma_{yy} = \displaystyle \sigma_{zz} = \displaystyle \tau_{yz} = 0$$

Hooke's Law for Three-Dimensional Isotropic Elastic Material
For the purposes of ascertaining strain from given stresses and vice versa, a 6x6 transformation matrix must be developed. But first, individual stress-strain relationships must be determined. From Hooke's Law the normal strains are:


 * $$ \epsilon_{xx} = \frac{\sigma_{xx}}{E} - \frac{\nu \sigma_{yy}}{E} - \frac{\nu \sigma_{zz}}{E} $$
 * $$ \epsilon_{yy} = \frac{\sigma_{yy}}{E} - \frac{\nu \sigma_{xx}}{E} - \frac{\nu \sigma_{zz}}{E} $$
 * $$ \epsilon_{zz} = \frac{\sigma_{zz}}{E} - \frac{\nu \sigma_{yy}}{E} - \frac{\nu \sigma_{xx}}{E} $$

while the shear strains are:


 * $$ \gamma_{xy} = 2\epsilon_{xy} = \frac{\tau_{xy}}{G} $$
 * $$ \gamma_{xz} = 2\epsilon_{xz} = \frac{\tau_{xz}}{G} $$
 * $$ \gamma_{yz} = 2\epsilon_{yz} = \frac{\tau_{yz}}{G} $$

The stresses and strains can be arranged in a column matrix as follows:


 * $$ (\epsilon_{ij})_{6x1} = \begin{bmatrix} \epsilon_{11} \\ \epsilon_{22} \\ \epsilon_{33} \\ \epsilon_{23} \\ \epsilon_{13} \\ \epsilon_{12} \end{bmatrix} $$  and   $$ (\sigma_{ij})_{6x1} = \begin{bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{33} \\ \sigma_{23} \\ \sigma_{13} \\ \sigma_{12} \end{bmatrix} $$

Combining the normal and shear strains into a transformation matrix yields the relations in the following HW problem solution:

MATLAB Airfoil Code Validation
Part of the first version of the MATLAB problem was to verify that the code was working properly by testing it both against a known area and also for replication of results. In the collapsible table below there is a "validation code" which calls part of the HW3 MATLAB code as a subfunction.

The resulting output from running the validation code is provided below.

The first section of the validation involved a "circular" airfoil. The first task was to determine the area of a circle of radius 1 meter at various numbers of contour segments. A table summarizing the above results is as follows:


 * {| border=1


 * bgcolor="red"|Number of Segments
 * bgcolor="#99ccff" style="text-align: center;" |20
 * bgcolor="#99ccff" style="text-align: center;" |40
 * bgcolor="#99ccff" style="text-align: center;" |60
 * bgcolor="#99ccff" style="text-align: center;" |80
 * bgcolor="#99ccff" style="text-align: center;" |100
 * bgcolor="red"|Calculated Area
 * bgcolor="#99ccff" |3.1045
 * bgcolor="#99ccff" |3.1285
 * bgcolor="#99ccff" |3.1344
 * bgcolor="#99ccff" |3.1369
 * bgcolor="#99ccff" |3.1383
 * }
 * }

The second task was to determine the $$ns_{max}\;$$, the number of segments at which the calculated area is within 1% of the actual area, which in this case is equal to 3.1416 square meters. This was done by calculating the area at various values of n and then comparing them. A plot of the results is given in Figure 9, where $$ns_{max}\;$$ and its corresponding area are denoted by a black square. For reference, the calculated value for $$ns_{max}\;$$ was found to be 23 segments, and the calculated area at this point was 3.1115 square meters, which differs from the actual value by .96%.

The final task for this section was to calculate the circle area at various locations of the origin point $$P_{0}\;$$. The three locations used were the coordinate system origin $$[0,0]\;$$, the circle center $$[1,0]\;$$, and an arbitrary external point $$[1,-4]\;$$. A table summarizing the results is shown below. As expected, the area is identical in all three cases.


 * {| border=1


 * bgcolor="red"|Coordinate of Origin Point
 * bgcolor="#99ccff" style="text-align: center;" |[0,0]
 * bgcolor="#99ccff" style="text-align: center;" |[1,0]
 * bgcolor="#99ccff" style="text-align: center;" |[1,-4]
 * bgcolor="red"|Calculated Area (square meters)
 * bgcolor="#99ccff" |3.1115
 * bgcolor="#99ccff" |3.1115
 * bgcolor="#99ccff" |3.1115
 * }
 * }

The second section of the validation involved a NACA 2415 airfoil of chord length c = 1 meter. The first task was to find $$ns_{max}\;$$ as with the circular airfoil. However, instead of comparing against a known area, the calculated area at a given number of segments is compared with the calculated area using one less segment, such that $$ns_{max}\;$$ is located when the difference between two successive area calculations is less than 1% of the current area calculation. A plot of the results is given in Figure 10. For reference, the calculated value for $$ns_{max}\;$$ was found to be 9 segments, and the calculated area at this point was 0.0244 square meters.

The second task was to show that the area calculation is valid at any choice of the origin point $$P_{0}\;$$. The three locations used were the coordinate system origin $$[0,0]\;$$, the airfoil centroid $$[.210,.008]\;$$ (previously calculated in HW3), and an arbitrary external point $$[.25,-1]\;$$. A table summarizing the results is shown below. Unlike with the circular case, the calculated areas are not identical for each case. However, a quick relative error calculation shows that the calculated values are all within 1% of each other, and such a small difference can be easily explained by the small numerical errors inherent in this kind of analysis.


 * {| border=1


 * bgcolor="red"|Coordinate of Origin Point
 * bgcolor="#99ccff" style="text-align: center;" |[0,0]
 * bgcolor="#99ccff" style="text-align: center;" |[.210,.008]
 * bgcolor="#99ccff" style="text-align: center;" |[.25,-1]
 * bgcolor="red"|Calculated Area (square meters)
 * bgcolor="#99ccff" style="text-align: center;" |.0244
 * bgcolor="#99ccff" style="text-align: center;" |.0246
 * bgcolor="#99ccff" style="text-align: center;" |.0247
 * }
 * }

Description of Airfoil Geometry and Equations


Figure 11 shows the geometry of a generic NACA 4-digit airfoil. The various quantities are defined as follows:


 * $$LE\;$$ - The point corresponding to the leading edge of the airfoil.
 * $$TE\;$$ - The point corresponding to the trailing edge of the airfoil.
 * $$c\;$$ - The chord length, describing the length of the straight line connecting the leading edge to the trailing edge.
 * $$m\;$$ - The maximum camber, corresponding to the first digit of the 4-digit series, and given in percentage of chord.
 * $$p\;$$ - The location of maximum camber, corresponding to the second digit of the 4-digit series, and given in tenths of chord.
 * $$t_{a}\;$$ - The maximum thickness, corresponding to the third and fourth digits of the 4-digit series, and given in percentage of chord.
 * $$y_{u}\;$$ - The contour defining the airfoil's upper surface.
 * $$y_{l}\;$$ - The contour defining the airfoil's lower surface.
 * $$MCL\;$$ - The contour defining the airfoil's mean camber line. Referred to as $$z_{c}\;$$ in the airfoil equations.

A set of general equations used to generate the various contours for an arbitrary NACA 4-digit airfoil can be found here.

Torsional Constant Comparison Between Single- and Multi-celled Airfoils
The new facet to the MATLAB assignment deals with a comparison of single and multi-celled airfoils, in terms of the torsional constant $$J\;$$. The code for the single-celled airfoil is in the table below:



The first portion of the assignment was to determine $$ns_{max}\;$$, as in previous assignments. The next collapsed table contains a "header" that can be applied to the above single celled airfoil analysis code to perform an iterative operation to generate an $$ns_{max}\;$$.The result was a $$ns_{max}\;$$ value of 12 segments and a $$J\;$$ value of $$4.7543x10^{-6} m^{4}\;$$. The second portion involved generating a graph of torsional constant vs. number of segments. The result is shown in Figure 12.

The second portion of the assignment was to determine $$J\;$$ for a NACA2415 airfoil with 3 cells. The MATLAB code used to carry out this task is given in a table below.



The plotted airfoil, with partitions at quarter and three-quarter chord, is shown in Figure 13. The calculated $$J\;$$ for this case, using $$ns = 1000\;$$ is $$5.2500x10^{-6} m^{4}\;$$.

Comparison:

The two calculated values for $$J\;$$ are very close to each other. The difference is likely a result of numerical errors during the generation of the torsional constant in the above MATLAB codes. Accounting for this effect, the overall result echoes the conclusion given in Chapter 3 of Sun: the torsional rigidity of a thin-walled, closed cross-section cannot be appreciably improved by compartmentalization.

MATLAB Code Certification
The following students are proficient in MATLAB programming:

Kevin Klauk Eas4200c.f08.vqcrew.E 22:01, 22 October 2008 (UTC)

Adam Edstrand Eas4200c.f08.vqcrew.a 05:09, 24 October 2008 (UTC)

Steven Hepsworth Eas4200c.f08.vqcrew.f 23:59, 23 October 2008 (UTC)

Darin Toscano Eas4200c.f08.vqcrew.d 07:28, 24 October 2008 (UTC)

Keith Javier Stober Eas4200c.f08.vqcrew.b 16:23, 24 October 2008 (UTC)

John Saxon Eas4200c.f08.vqcrew.c 20:41, 24 October 2008 (UTC)

Contributing Team Members
The following students contributed to this report:

Kevin Klauk Eas4200c.f08.vqcrew.E 22:01, 22 October 2008 (UTC)

Adam Edstrand Eas4200c.f08.vqcrew.a 05:09, 24 October 2008 (UTC)

Steven Hepsworth Eas4200c.f08.vqcrew.f 23:59, 23 October 2008 (UTC)

Darin Toscano Eas4200c.f08.vqcrew.d 07:29, 24 October 2008 (UTC)

Keith Javier Stober Eas4200c.f08.vqcrew.b 16:24, 24 October 2008 (UTC)

John Saxon Eas4200c.f08.vqcrew.c 20:42, 24 October 2008 (UTC)