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

See my comments below. Eml4500.f08 21:54, 10 December 2008 (UTC)

Non-Symmetric Thin Walled Cross-Section

Fig 1. - Non-symmetrical thin-walled cross section with stringers varying in area.

Now that the formulas have been laid out for non-symmetric case, as seen in HW6, it is now time to look at a problem utilizing these equations. The non-symmetric case can be seen in Fig 1. This case is vary similar to the symmetrical cross-section in example problem 5.2 in the book, however the areas will vary by stringer for the non-symmetrical case. The figure also shows the location of the centroid which was calculated in the HW section below, and the shear flow across each skin panel.

$A_{1}=A\;$

$A_{2}=2A\;$

$A_{3}=3A\;$

$A_{4}=4A\;$

In evaluating Fig 1., the mean value theorem will be used in analyzing the first moment of the area, Q_y and Q_z, so the equations can be modified to the following:

$Q_{z}=\int _{A}{\bar {y}}\,dA={\bar {y}}A$

$Q_{y}=\int _{A}{\bar {z}}\,dA={\bar {z}}A$

The Area to be used for the overall moments will have to be the sum of the area of the components. Neglecting the skin and spar webs the equation for the overall area the sum of the area of $A_{i}\;$

$A=\sum _{i=1}^{4}{A_{i}}\;$

From the first moment of area equations and the known area it is now possible to examine the shear flow along the path s. Looking at Fig 1., V_z is zero for this case so the shear flow equations becomes the following:

$q(s)=(k_{yz}Q_{z}^{(s)}-k_{z}Q_{y}^{(s)})V_{z}$

Now if we break down the shear flow equation and examine each variable it can be determined the each variable is actually independent of the past s.

$\left.{\begin{matrix}1)&V_{z}&\\2)&k_{yz}&k_{z}\\3)&Q_{z}&Q_{y}\end{matrix}}\right\}$ independent of s

From the two first moment of inertia equations it can also be seen the they are constant between two stringers. This is due to all areas are concentrated of stringers, thus constant and independent of s between the stringers.

These relations yield the final assumption that the shear flow q(s) is constant between two stringers, but q(s) would increment (jump) when crossing over a stringer.

Now that we have the following relation with the stringers, it is now time to lay out a game plan to solve for the shear flows along the path s.

Step 1: Find $({\bar {y}},{\bar {z}})$

Step 2: Find $I_{y},I_{z},I_{yz}\;$

Step 3: Find $k_{y},k_{z},andk_{y}z\;$

Step 4: Follow path "s" to find $q_{1-2},q_{2-3},q_{3-4}\;$

$q_{1-2}\;$ represents the shear flow in skin panel 1-2. This can be represented in the same form as the shear flow along any path.

$q_{12}=(k_{yz}Q_{z}^{12}-k_{z}Q_{y}^{12})V_{z}$

Where the first moment area of inertia is local around each panel. Where $y_{1}$ is the y component of the string with respect to the origin at the centroid, and similarly for $z_{1}$

$Q_{z}^{12}=y_{1}A_{1}$

$Q_{y}^{12}=z_{1}A_{1}$

Fig 2. - Non-symmetrical thin-walled cross section showing various areas and paths chosen.

A quick look for reference of the following skin panel 2-3 shows what the first moment of inertia would resemble in that case.

$Q_{z}^{23}=y_{1}A_{1}+y_{2}A_{2}$

A look at a different type of non-symmetrical case might help better explain different approaches that could be taken when looking for the shear flow. Fig 2. shows different paths and areas for each path that could be used for solving such a shape.

The area for this case that was split into two section is the sum of the two sections:

$A=A_{\hat {s}}+A_{{\hat {s}},L}$

HW Problem - Find the centroid for Fig. 1

For this problem the areas stated above were used in the calculations of the centroid. $A_{1}=A\;$ $A_{2}=2A\;$ $A_{3}=3A\;$ $A_{4}=4A\;$ The following equations were found in Mechanic of Aircraft Structures , for calculating the location of the centroid. ${\bar {y_{c}}}={\frac {\sum {{\bar {y_{i}}}A_{i}}}{\sum {A_{i}}}}$ ${\bar {z_{c}}}={\frac {\sum {{\bar {z_{i}}}A_{i}}}{\sum {A_{i}}}}$ Now to find $({\bar {y_{c}}},{\bar {z_{c}}})$. For calculation purposes the length of each skin panel will be taken to be 1m and $A=10^{-}2m^{2}$. ${\bar {y_{c}}}={\frac {\sum {{\bar {y_{i}}}A_{i}}}{\sum {A_{i}}}}={\frac {1(A_{1})+1(A_{4})}{A_{1}+A_{2}+A_{3}+A_{4}}}=0.5m$ ${\bar {y_{c}}}={\frac {\sum {{\bar {y_{i}}}A_{i}}}{\sum {A_{i}}}}={\frac {1(A_{1})+1(A_{2})}{A_{1}+A_{2}+A_{3}+A_{4}}}=0.3m$ Thus the centroid is located 0.5 meters to the right of the origin, or half-way between the stringers, and 0.3 m above the origin. $({\bar {y_{c}}},{\bar {z_{c}}})=(0.5,0.3)$

Mini Plan

Single-Cell Sections

- without stringers

- with stringers

Multi-Cell Sections

- without stringers

- with stringers

Single-Cell Section Without Stringers

Without Stringers

Fig.3 - A thin-walled cross-section without stringers being acted upon by a shear flow.

Question: Can the setup shown in Figure 3 resist $\displaystyle V_{z}$?

Answer: No.

$R^{z}=R_{AB}^{z}+R_{BA}^{z}$

$\displaystyle R^{z}$ is the total resultant in the z-direction

$R_{AB}^{z}$ is the reslutant of q in AB

$R_{BA}^{z}$ is the resultant of q in BA

$R_{AB}^{z}=-q{\bar {A^{'}B^{'}}}=R_{BA}^{z}$

so, $\displaystyle R^{z}=0$

Therefore, the configuration can not resist any $\displaystyle V_{z}$.

With Stringers

Fig. 4 - A single-cell section with stringers being acted upon by different shear flows in different sections.

Fig. 5 - P1: A single section without stringers to allow for the shear flow in the case with stringers.

Fig. 6 - P2: A single cross-section with stringers and a wall cut allowing for the determination of shear flow.

Figure 4 shows that q(s) is piecewise constant with respect to the curvilinear coordinate, s. (i.e. $\displaystyle q_{ij}$ = constant within each panel)

However, the value for q is not the same throughout all of the panels.

$q_{12}\neq q_{23}\neq q_{31}$

The presence of the stringers causes a non-constant q.

$R^{z}=V_{z}\neq 0$

Superposition can be used to solve for the values of $\displaystyle q_{ij}$ in each section by solving for values of q in the case without stringers, Figure 5, and for the case including an imaginary cut, Figure 6, in one of the sections. This is a valid use of superposition due to the linearity of the relationship.

$q_{12}=q+{\tilde {q}}_{12}$

$q_{23}=q+{\tilde {q}}_{23}$

$q_{31}=q+{\tilde {q}}_{31}$

Analysis Algorithm

Objectives: One unknown q ($\displaystyle {\tilde {q_{ij}}}$ are known after solving P1), need to solve one equation for one unknown.

Method: Knowing the values for $\displaystyle V_{y}$ and $\displaystyle V_{z}$

1) Solve P2 for $\displaystyle {\tilde {q_{12}}}$ and ${\tilde {q_{23}}}$

2) Moment Equation - take the moment about any point in the plane (y,z)

a) Superposition: $\displaystyle q_{ij}=q+{\tilde {q_{ij}}}$

b) Select point $\displaystyle {\bar {O}}$ in the plane (y,z)

Multi-Cell Sections

Without Stringers

Fig. 7 - A multi-cell section without stingers being acted upon by several shear flows.

The same type of analysis can be performed for multi-cell sections as for single cell sections. Figure 7 shows the multi-cell case without stringers. The resultant of all of the shear flows in each cell can be calculated by the sum of the resultants in every cell.

$\displaystyle R^{z}=\sum _{i=1}^{N_{cells}}R_{i}^{z}$

Similar to the single case, the resultant for each cell can be found, therefore yielding the result for the entire multi-cell section. Once again, this value is equal to zero. Without the support of the stringers, this setup can not resist $\displaystyle V_{z}$.

$\displaystyle R_{i}^{z}=0\Rightarrow R^{z}=0$

With Stringers

Fig. 8 - A multi-cell section with stringers.

Fig. 9 - P1: A multi-cell section without stringers to aid in the determination of shear flow in the case with stringers

Fig. 10 - P2: A multi-cell section with stingers and wall cuts to help solve for the case with stringer.

The true shear flow, $\displaystyle q_{ij}$, is equal to the sum of the closed cell constant shear flow, q, and the open cell piecewise constant shear flow, $\displaystyle {\tilde {q_{ij}}}$.

$\displaystyle q_{ij}=q+{\tilde {q_{ij}}}$

In order to solve for all of the true shear flows, imaginary cuts must be made into some of the walls. This must be done carefully to ensure that no stringers are left unattached to the rest of the assembly. Figure 8 shows the overall problem including no cuts and all of the stringers. Figure 9 is displaying the case in which the stringers are removed so that the closed cell shear flow can be calculated. Figure 10 represents the case in which the cuts have been made and the open cell piecewise constant shear flows can be found. The total resultant of the multi-cell case with stringers can be calculated by the sum of the resultants of P1 and P2.

$\displaystyle R^{z}=R^{z_{1}}+R^{z_{2}}$

Recall: $\displaystyle R^{z_{1}}=0$

Therefore: $\displaystyle R^{z}=R^{z_{2}}=V_{Z}\neq 0$

HW Problem - Derive further required Q values

We have already solved for the values of $Q_{z}^{2}\,$ and $Q_{z}^{3}\,$ but what about the values of their counterparts, $Q_{y}^{2}\,$ and $Q_{y}^{3}\,$? It turns out we will use the very same formulas to compute these! That is equations 5.5 and 5.6, which tell us that: $q_{s}=-(k_{y}V_{y}-k_{yz}V_{z})Q_{z}-(k_{z}V_{z}-k_{yz}V_{y})Q_{y}\,$ $Q_{z}=\int _{As}\int ydA\,$ $Q_{y}=\int _{As}\int zdA\,$ In much the same way as before, we are able to utilize the final segment of these equations to derive Q values for both the second and third stringers... $Q_{y}^{2}=z_{2}A_{2}\,$ $Q_{y}^{3}=z_{3}A_{3}\,$

Solving P2 - Euler Cut Principle

The Euler Cut Principle involves finding the equilibrium of each stringer. This is basically a statics problem where the sum of the forces, which are given through the shear flows, are equal to zero. Stringer 3 is shown in Figure 11 in two dimensions and again in three dimensions in Figure 12, with the stringer expanded into the x-direction. The sum of the forces in the x-direction can be expressed as follows:

Fig 11. 2-D drawing of Stringer 3

Fig 12. 3-D drawing of Stringer 3

Fig 13. General form as consequence of first method

Fig 14. Stringer j with left hand side in opposite direction

Fig 15. Airfoil with cuts to isolate Stringer 3

Fig 16. 2-D view of Stringer 3 with cuts

$\sum {F_{x}}=0=\int _{A_{3}}^{}{\left[\sigma _{xx}(x+dx)-\sigma _{xx}(x)\right]}dA_{3}=\int _{A_{3}}^{}{\left[{\frac {d\sigma _{xx}}{dx}}dx+HOT\right]}dA_{3}=\left[{\tilde {q_{31}}}-{\tilde {q_{23}}}-{\tilde {q_{43}}}\right]dx$

Thus we solve for ${\tilde {q_{31}}}$ as:

${\tilde {q_{31}}}={\tilde {q_{23}}}+{\tilde {q_{43}}}+q^{(3)}$

Where $q^{(3)}\;$ is the contribution to shear flow by Stringer 3:

$q^{(3)}=-\int _{A_{3}}^{}{\frac {d\sigma _{xx}}{dx}}dA_{3}$

The integral can be solved for $q^{(3)}\;$ as:

$q^{(3)}=-\left(k_{y}V_{y}-k_{yz}V_{z}\right)Q_{z}^{(3)}-\left(k_{z}V_{z}-k_{yz}V_{y}\right)Q_{y}^{(3)}$

where:

$V_{y}={\frac {dM_{z}}{dx}}$

$V_{z}={\frac {dM_{y}}{dx}}$

$Q_{z}^{(3)}=\int _{A_{3}}^{}{ydA_{3}}$

$Q_{y}^{(3)}=\int _{A_{3}}^{}{zdA_{3}}$

The same procedure can be used to solve for Stringer 2:

${\tilde {q_{23}}}=\underbrace {\tilde {q_{12}}} _{flowin}-\underbrace {\tilde {q_{24}}} _{flowout}+q^{(2)}$

The value of $q^{(2)}\;$ can be computed as done with Stringer 3 with the values of $Q_{z}^{(2)}$ and $Q_{y}^{(2)}$ as:

$Q_{z}^{(2)}=y_{2}A_{2}$

$Q_{y}^{(2)}=z_{2}A_{2}$

where $y_{2}\;$ and $z_{2}\;$ are the y and z coordinates of Stringer 2, and $A_{2}\;$ is the area of Stringer 2.

The same procedure can be followed on Stringer 4 to find $q^{(4)}\;$ and ${\tilde {q_{24}}}$. Hence, we can also deduce ${\tilde {q_{31}}}$.

Using the principles of superposition, the values of q can be determined.

$q_{ij}={\tilde {q}}_{ij}+q_{k}$

On the right side of the equations, the tilde terms are the knowns and the non-tilde terms are unknowns.

$q_{12}={\tilde {q}}_{12}+q_{1}$

$q_{23}={\tilde {q}}_{23}+q_{1}$

$q_{31}={\tilde {q}}_{31}+q_{1}$

$q_{24}={\tilde {q}}_{24}+q_{2}$

$q_{43}={\tilde {q}}_{43}+q_{2}-q_{3}$

$q_{41}={\tilde {q}}_{41}+q_{3}$

These equations reveal that there are three unknowns, $q_{1},q_{2},q_{3}\;$, therefore three equations are needed to solve for these values to complete the problem. The three equations to be use are:
1) Moment equation: Take the moment of $V_{y}\;$, $V_{z}\;$, and $\left\{q_{12},...,q_{41}\right\}$ about any convenient point (usually where the lines of action of $V_{y}\;$ and $V_{z}\;$ intersect).
2) $\theta _{1}=\theta _{2}\;$
3) $\theta _{2}=\theta _{3}\;$
Equations 2 and 3 are the compatibility equations.

To recap, in order to solve P2 for each cell you must follow path "$s_{i}\;$" and then determine the equilibrium of each stringer on the path. There are two ways to accomplish this:
1) Complete method using FBD as shown above
2) A consequence of the first method, as exemplified in Figure 13:

${\tilde {q_{j6}}}={\tilde {q_{j2}}}-{\tilde {q_{j5}}}-{\tilde {q_{j8}}}+q^{(j)}$

Also, the left side of the equation simply becomes negative if the direction of ${\tilde {q_{j6}}}\;$ changes (and becomes ${\tilde {q_{6j}}}$), as shown in Figure 14.

${\tilde {q_{6j}}}={\tilde {q_{j2}}}-{\tilde {q_{j5}}}-{\tilde {q_{j8}}}+q^{(j)}$

Using the Euler Cut Principle brings forth a question, what happens if the object is cut such that one stringer is isolated as shown in Figure 15? Due to the cuts made, ${\tilde {q_{23}}}={\tilde {q_{31}}}={\tilde {q_{34}}}=0$. We then sum the forces acting on the FBD shown in Figure 16.

$\underbrace {\tilde {q_{31}}} _{0}=\underbrace {\tilde {q_{23}}} _{0}-\underbrace {\tilde {q_{34}}} _{0}+q^{(3)}$

This shows that $q^{(3)}=0\;$, which is not possible. Therefore, when using the Euler Cut method a stringer must not be isolated from the rest of the body.

Shear Buckling

Expressing ${\begin{Bmatrix}C_{22},C_{13},C_{31},C_{33}\end{Bmatrix}}$ in terms of $C_{11}$ and for $\vartheta =1.5$ we can accomplish:

1) Finding $\lambda$ for $\vartheta =1.5$

2) Evaluate numerically for $k_{sys}$ in equation 30

3) Find $[k_{ij}]$

${\begin{bmatrix}k_{12}&k_{23}&\cdots &k_{25}\\\vdots &\ddots &&\vdots \\\vdots &&\ddots &\vdots \\k_{32}&\cdots &\cdots &k_{55}\end{bmatrix}}{\begin{bmatrix}C_{22}\\\vdots \\\vdots \\C_{33}\end{bmatrix}}={\begin{bmatrix}-{\frac {u_{z}}{a}}C_{11}\\0\\0\\0\end{bmatrix}}$

Then solve ${\begin{Bmatrix}C_{22},C_{13},C_{31},C_{33}\end{Bmatrix}}$ intersect with $C_{11}$. If the $K$ matrix is expressed as ${\underline {\overline {K}}}$ then the inverse of said matrix is ${\underline {\overline {K}}}^{-1}$. The values for ${\begin{Bmatrix}C_{22},C_{13},C_{31},C_{33}\end{Bmatrix}}$ can be related as

${\begin{Bmatrix}C_{22},C_{13},C_{31},C_{33}\end{Bmatrix}}={\underline {\overline {K}}}^{-1}\ {\begin{bmatrix}-{\frac {u_{z}}{a}}C_{11}\\0\\0\\0\end{bmatrix}}$

Resulting in

$u_{z}=C_{11}sin({\frac {\pi x}{a}})sin({\frac {\pi y}{b}})+C_{22}sin({\frac {2\pi x}{a}})sin({\frac {2\pi y}{b}})+C_{13}sin({\frac {\pi x}{a}})sin({\frac {3\pi y}{b}})+C_{31}sin({\frac {3\pi x}{a}})sin({\frac {\pi y}{b}})+C_{33}sin({\frac {3\pi x}{a}})sin({\frac {3\pi y}{b}})$

HW Problem - Finding $q_{43}$ and $q_{41}$
In class it was illustrated that in the asymmetric thin-walled cross-section after the Euler cuts were made that $q_{24}$ could be related as $q_{24}={\tilde {q}}_{24}+q_{2}$ It was further required that the values of $q_{43}$ and $q_{41}$ were needed to be found as well. It can be shown that $q_{43}={\tilde {q}}_{43}+q_{4}$ $q_{41}={\tilde {q}}_{41}+q_{4}$

Answer to Jared's Question of Airfoil Cell Lacking Exposed Surface Area - Part 2

In the case of an airfoil lacking exposed surface area, the resulting equation of shear flow would be as follows:

$0=q^{(1)}+q^{(2)}+q^{(4)}$

Which cannot be true or even possible, since:

$0=\sum _{e=1}^{4}q^{(e)}$

This equation can also be expressed as:

$0=q^{(1)}+q^{(2)}+q^{(3)}+q^{(4)}\;\therefore \;q^{(1)}+q^{(2)}+q^{(4)}=-q^{(3)}\neq 0$

To prove that the sum of all the shear flows is zero (as shown above) the following equation, derived earlier, can be utilized:

$q^{(e)}=n_{z}Q_{z}^{(e)}+n_{y}Q_{y}^{(e)}$

It should be noted that n in the above equation does not represent a unit normal directional component, as was the case previously. In the above equation:

$n_{z}=-(k_{y}V_{y}-k_{yz}V_{z})$

$n_{y}=-(k_{z}V_{z}-k_{yz}V_{z})$

Now, the above expression for shear flow, q, can be used to reinforce the fact that the sum of all the individual shear flows is zero.

$\sum _{e=1}^{4}q^{(e)}=n_{z}\sum _{e=1}^{4}Q_{z}^{(e)}+n_{y}\sum _{e=1}^{4}Q_{y}^{(e)}$

Using the definition of $Q_{z}$ and $Q_{y}$ along with the Mean Value Theorem, it can be shown that:

$\sum _{e=1}^{4}Q_{z}^{(e)}=0\;and\;\sum _{e=1}^{4}Q_{y}^{(e)}=0\;\therefore \;\sum _{e=1}^{4}q^{(e)}=0$

Read and Report: Section 4.2 Bidirectional Bending

When considering a beam characterized by bidirectional bending, the approximate displacement expansions, u, v and w, are as follows:

$u=u_{0}(x)+z\psi _{y}(x)+y\psi _{z}(x)$

$v=v_{0}(x)$

$w=w_{0}(x)$

In the above equations, $\psi _{y}$ and $\psi _{z}$ represent rotations about the y and z axes, respectively, using the right-hand rule about the respective axes to determine the sign of each. Consequently, the strains are as follows:

$\epsilon _{xx}={\frac {\delta u}{\delta x}}={\frac {du_{0}}{dx}}+z{\frac {d\psi _{y}}{dx}}+y{\frac {d\psi _{z}}{dx}}$

$\gamma _{xy}={\frac {\delta v}{\delta x}}+{\frac {\delta u}{\delta y}}={\frac {dv_{0}}{dx}}+\psi _{z}$

$\gamma _{xz}={\frac {\delta w}{\delta x}}+{\frac {\delta u}{\delta z}}={\frac {dw_{0}}{dx}}+\psi _{y}$

Using the assumption that $\gamma _{xy}=0$ and $\gamma _{xz}=0$ simplifies the last two expressions to:

$\psi _{z}=-{\frac {dv_{0}}{dx}}$

$\psi _{y}=-{\frac {dw_{0}}{dx}}$

which when subsituted into the equation for strain in the x direction yield:

$\epsilon _{xx}={\frac {du_{0}}{dx}}-y{\frac {d^{2}v_{0}}{dx^{2}}}-z{\frac {d^{2}w_{0}}{dx^{2}}}$

If $N_{x}=0$, then ${\frac {du_{0}}{dx}}=0$, further simplifying the equation of x-axis strain to:

$\epsilon _{xx}=-y{\frac {d^{2}v_{0}}{dx^{2}}}-z{\frac {d^{2}w_{0}}{dx^{2}}}$

It is next desirable to apply this result to the y and z bending moment equations. The definition of these moments is as follows:

$M_{y}=\int _{A}\int z\sigma _{xx}dA=-E\int _{A}\int (yz{\frac {d^{2}v_{0}}{dx^{2}}}+z^{2}{\frac {d^{2}w_{0}}{dx^{2}}})dA=-EI_{yz}{\frac {d^{2}v_{0}}{dx^{2}}}-EI_{y}{\frac {d^{2}w_{0}}{dx^{2}}}$

$M_{z}=\int _{A}\int y\sigma _{xx}dA=-EI_{yz}{\frac {d^{2}w_{0}}{dx^{2}}}-EI_{z}{\frac {d^{2}v_{0}}{dx^{2}}}$

where

$I_{y}=\int _{A}\int z^{2}dA$

$I_{z}=\int _{A}\int y^{2}dA$

$I_{yz}=\int _{A}\int yzdA$

In the definition of bending moment equations, the quadratic terms can be expressed as the curvature, i.e.:

$\chi _{y}={\frac {d^{2}v_{0}}{dx^{2}}}\;and\;\chi _{z}={\frac {d^{2}w_{0}}{dx^{2}}}$

Now, the equations can be rewritten in matrix form:

${\begin{bmatrix}M_{y}\\M_{z}\end{bmatrix}}={\begin{bmatrix}I_{y}&I_{yz}\\I_{yz}&I_{z}\end{bmatrix}}{\begin{bmatrix}-E\chi _{z}\\-E\chi _{y}\end{bmatrix}}=\mathbf {I} {\begin{bmatrix}-E\chi _{z}\\-E\chi _{y}\end{bmatrix}}$

$\epsilon _{xx}=-y\chi _{y}-z\chi _{z}={\begin{bmatrix}z&y\end{bmatrix}}{\begin{bmatrix}-\chi _{z}\\-\chi _{y}\end{bmatrix}}$

$\sigma _{xx}=E\epsilon _{xx}=E{\begin{bmatrix}z&y\end{bmatrix}}{\begin{bmatrix}-\chi _{z}\\-\chi _{y}\end{bmatrix}}=E{\begin{bmatrix}z&y\end{bmatrix}}{\frac {1}{E}}\mathbf {I^{-1}} {\begin{bmatrix}M_{y}\\M_{z}\end{bmatrix}}={\begin{bmatrix}z&y\end{bmatrix}}\mathbf {I^{-1}} {\begin{bmatrix}M_{y}\\M_{z}\end{bmatrix}}$

where

$\mathbf {I^{-1}} ={\frac {1}{D}}{\begin{bmatrix}I_{z}&-I_{yz}\\-I_{yz}&I_{y}\end{bmatrix}}$

and

$D=I_{y}I_{z}-(I_{yz})^{2}$

Similarly, the shear flow equation can be converted to matrix form:

$q=-\int _{A}\int {\frac {d\sigma _{xx}}{dx}}dA$

where

${\frac {d\sigma _{xx}}{dx}}={\begin{bmatrix}z&y\end{bmatrix}}\mathbf {I^{-1}} {\begin{bmatrix}{\frac {dM_{y}}{dx}}\\{\frac {dM_{z}}{dx}}\end{bmatrix}}$

Recall that:

$V_{y}={\frac {dM_{y}}{dx}}$

$V_{z}={\frac {dM_{z}}{dx}}$

and

$Q_{y}=\int _{A}\int zdA$

$Q_{z}=\int _{A}\int ydA$

MATLAB Coding Problem Continued

Full MATLAB Code - Single Celled Airfoil, Shear Analysis

function vqcrew_shear_singlecell(c,m,p,t,ns) close all % Define given parameters c = c+.001; ts = .002; %skin thickness tw = .003; %spar thickness tstr = .005; %stringer thickness A_B = 2e-4; %stringer areas A_E = 2e-4; A_H = 1e-4; A_F = 1e-4; A_T = A_B+A_E+A_H+A_F; %total stringer area M_y = -1250; %y direction bending moment M_z = 500; %z direction bending moment V_y = 10000; %horizontal shear force V_z = 5000; %vertical shear force %Generate airfoil contour and contour length [yu,zu,yl,zl] = airfoil(c,m,p,t,ns); su = pathlength(yu,zu); sl = pathlength(yl,zl); s = su+sl; 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; 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; 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; % Define point B stringer endpoints and stringer centroid Lbrack = .5*A_B/tstr; ay1 = min(abs(yu-.25*c-.5*tw)); by1 = min(abs(yu-.25*c+.5*tw)); ay2 = min(abs(yu-.25*c-.5*tw-.5*Lbrack)); by2 = min(abs(yu-.25*c+.5*tw+.5*Lbrack)); j=1; while abs(yu(j)-.25*c-.5*tw)~=ay1 j = j+1; end jb1 = j; j=1; while abs(yu(j)-.25*c+.5*tw)~=by1 j = j+1; end jb2 = j; j=1; while abs(yu(j)-.25*c-.5*tw-.5*Lbrack)~=ay2 j = j+1; end jb3 = j; j=1; while abs(yu(j)-.25*c+.5*tw+.5*Lbrack)~=by2 j = j+1; end jb4 = j; Bpoint1a = [yu(jb2);zu(jb2)-.5*ts]; Bpoint2a = [yu(jb4);zu(jb4)-.5*ts]; Bpoint3a = [Bpoint1a(1);Bpoint1a(2)-.5*Lbrack]; Bpoint4a = [Bpoint3a(1)-tstr;Bpoint3a(2)]; Bpoint5a = [Bpoint2a(1);Bpoint2a(2)-tstr]; Bpoint6a = [Bpoint4a(1);Bpoint5a(2)]; Bpoint1b = [yu(jb1);zu(jb1)-.5*ts]; Bpoint2b = [yu(jb3);zu(jb3)-.5*ts]; Bpoint3b = [Bpoint1b(1);Bpoint1b(2)-.5*Lbrack]; Bpoint4b = [Bpoint3b(1)+tstr;Bpoint3b(2)]; Bpoint5b = [Bpoint2b(1);Bpoint2b(2)-tstr]; Bpoint6b = [Bpoint4b(1);Bpoint5b(2)]; Bcent(1,1) = (Bpoint2a(1) + Bpoint2b(1) + Bpoint5a(1) + Bpoint5b(1) + ... Bpoint6a(1) + Bpoint6b(1) + Bpoint4a(1) + Bpoint4b(1))/8; Bcent(2,1) = (Bpoint2a(2) + Bpoint2b(2) + Bpoint5a(2) + Bpoint5b(2) + ... Bpoint6a(2) + Bpoint6b(2) + Bpoint4a(2) + Bpoint4b(2))/8; % Define point E stringer endpoints and stringer centroid Lbrack = .5*A_E/tstr; ay1 = min(abs(yl-.25*c-.5*tw)); by1 = min(abs(yl-.25*c+.5*tw)); ay2 = min(abs(yl-.25*c-.5*tw-.5*Lbrack)); by2 = min(abs(yl-.25*c+.5*tw+.5*Lbrack)); j=1; while abs(yl(j)-.25*c-.5*tw)~=ay1 j = j+1; end je1 = j; j=1; while abs(yl(j)-.25*c+.5*tw)~=by1 j = j+1; end je2 = j; j=1; while abs(yl(j)-.25*c-.5*tw-.5*Lbrack)~=ay2 j = j+1; end je3 = j; j=1; while abs(yl(j)-.25*c+.5*tw+.5*Lbrack)~=by2 j = j+1; end je4 = j; Epoint1a = [yl(je2);zl(je2)+.5*ts]; Epoint2a = [yl(je4);zl(je4)+.5*ts]; Epoint3a = [Epoint1a(1);Epoint1a(2)+.5*Lbrack]; Epoint4a = [Epoint3a(1)-tstr;Epoint3a(2)]; Epoint5a = [Epoint2a(1);Epoint2a(2)+tstr]; Epoint6a = [Epoint4a(1);Epoint5a(2)]; Epoint1b = [yl(je1);zl(je1)+.5*ts]; Epoint2b = [yl(je3);zl(je3)+.5*ts]; Epoint3b = [Epoint1b(1);Epoint1b(2)+.5*Lbrack]; Epoint4b = [Epoint3b(1)+tstr;Epoint3b(2)]; Epoint5b = [Epoint2b(1);Epoint2b(2)+tstr]; Epoint6b = [Epoint4b(1);Epoint5b(2)]; Ecent(1,1) = (Epoint2a(1) + Epoint2b(1) + Epoint5a(1) + Epoint5b(1) + ... Epoint6a(1) + Epoint6b(1) + Epoint4a(1) + Epoint4b(1))/8; Ecent(2,1) = (Epoint2a(2) + Epoint2b(2) + Epoint5a(2) + Epoint5b(2) + ... Epoint6a(2) + Epoint6b(2) + Epoint4a(2) + Epoint4b(2))/8; % Define point H stringer endpoints and stringer centroid Lbrack = .5*A_H/tstr; ay1 = min(abs(yl-.75*c-.5*tw)); by1 = min(abs(yl-.75*c+.5*tw)); ay2 = min(abs(yl-.75*c-.5*tw-.5*Lbrack)); by2 = min(abs(yl-.75*c+.5*tw+.5*Lbrack)); j=1; while abs(yl(j)-.75*c-.5*tw)~=ay1 j = j+1; end jh1 = j; j=1; while abs(yl(j)-.75*c+.5*tw)~=by1 j = j+1; end jh2 = j; j=1; while abs(yl(j)-.75*c-.5*tw-.5*Lbrack)~=ay2 j = j+1; end jh3 = j; j=1; while abs(yl(j)-.75*c+.5*tw+.5*Lbrack)~=by2 j = j+1; end jh4 = j; Hpoint1a = [yl(jh2);zl(jh2)+.5*ts]; Hpoint2a = [yl(jh4);zl(jh4)+.5*ts]; Hpoint3a = [Hpoint1a(1);Hpoint1a(2)+tstr]; Hpoint4a = [Hpoint3a(1)-tstr;Hpoint3a(2)]; Hpoint1b = [yl(jh1);zl(jh1)+.5*ts]; Hpoint2b = [yl(jh3);zl(jh3)+.5*ts]; Hpoint3b = [Hpoint1b(1);Hpoint1b(2)+tstr]; Hpoint4b = [Hpoint3b(1)+tstr;Hpoint3b(2)]; Hcent(1,1) = (Hpoint2a(1) + Hpoint2b(1) + Hpoint4a(1) + Hpoint4b(1))/4; Hcent(2,1) = (Hpoint2a(2) + Hpoint2b(2) + Hpoint4a(2) + Hpoint4b(2))/4; % Define point F stringer endpoints and stringer centroid Lbrack = .5*A_F/tstr; ay1 = min(abs(yu-.75*c-.5*tw)); by1 = min(abs(yu-.75*c+.5*tw)); ay2 = min(abs(yu-.75*c-.5*tw-.5*Lbrack)); by2 = min(abs(yu-.75*c+.5*tw+.5*Lbrack)); j=1; while abs(yu(j)-.75*c-.5*tw)~=ay1 j = j+1; end jf1 = j; j=1; while abs(yu(j)-.75*c+.5*tw)~=by1 j = j+1; end jf2 = j; j=1; while abs(yu(j)-.75*c-.5*tw-.5*Lbrack)~=ay2 j = j+1; end jf3 = j; j=1; while abs(yu(j)-.75*c+.5*tw+.5*Lbrack)~=by2 j = j+1; end jf4 = j; Fpoint1a = [yu(jf2);zu(jf2)-.5*ts]; Fpoint2a = [yu(jf4);zu(jf4)-.5*ts]; Fpoint3a = [Fpoint1a(1);Fpoint1a(2)-tstr]; Fpoint4a = [Fpoint3a(1)-tstr;Fpoint3a(2)]; Fpoint1b = [yu(jf1);zu(jf1)-.5*ts]; Fpoint2b = [yu(jf3);zu(jf3)-.5*ts]; Fpoint3b = [Fpoint1b(1);Fpoint1b(2)-tstr]; Fpoint4b = [Fpoint3b(1)+tstr;Fpoint3b(2)]; Fcent(1,1) = (Fpoint2a(1) + Fpoint2b(1) + Fpoint4a(1) + Fpoint4b(1))/4; Fcent(2,1) = (Fpoint2a(2) + Fpoint2b(2) + Fpoint4a(2) + Fpoint4b(2))/4; %Define point D pointD(1,1) = (Bcent(1,1) + Ecent(1,1))/2; pointD(2,1) = (Bcent(2,1) + Ecent(2,1))/2; %Find overall stringer centroid s_cent(1,1) = (Fcent(1,1)*A_F + Bcent(1,1)*A_B + Ecent(1,1)*A_E + ... Hcent(1,1)*A_H)/A_T; s_cent(2,1) = (Fcent(2,1)*A_F + Bcent(2,1)*A_B + Ecent(2,1)*A_E + ... Hcent(2,1)*A_H)/A_T; %Find Skin Area and Centroid A_skin = s*ts; sk_cent = centroid(yu,zu,yl,zl,ns); % Calculate Wing Centroid A_tot = A_T + A_skin; % Total wing component area cw(1,1) = (s_cent(1)*A_T + sk_cent(1)*A_skin)/(A_tot); cw(2,1) = (s_cent(2)*A_T + sk_cent(2)*A_skin)/(A_tot); I_22 = (A_B*Bcent(2)^2 + A_B*(Bcent(1)-cw(1))^2) + ... (A_E*Ecent(2)^2 + A_E*(Ecent(1)-cw(1))^2) + ... (A_H*Hcent(2)^2 + A_H*(Hcent(1)-cw(1))^2) + ... (A_F*Fcent(2)^2 + A_F*(Fcent(1)-cw(1))^2) + ... (A_skin*sk_cent(2)^2 + A_skin*(sk_cent(1)-cw(1))^2); I_33 = (A_B*Bcent(1)^2 + A_B*(Bcent(2)-cw(2))^2) + ... (A_E*Ecent(1)^2 + A_E*(Ecent(2)-cw(2))^2) + ... (A_H*Hcent(1)^2 + A_H*(Hcent(2)-cw(2))^2) + ... (A_F*Fcent(1)^2 + A_F*(Fcent(2)-cw(2))^2) + ... (A_skin*sk_cent(1)^2 + A_skin*(sk_cent(2)-cw(2))^2); I_23 = (A_B*(Bcent(2)-cw(2))*(Bcent(1)-cw(1))) + ... (A_E*(Ecent(2)-cw(2))*(Ecent(1)-cw(1))) + ... (A_H*(Hcent(2)-cw(2))*(Hcent(1)-cw(1))) + ... (A_F*(Fcent(2)-cw(2))*(Fcent(1)-cw(1))) + ... (A_skin*(sk_cent(2)-cw(2))*(sk_cent(1)-cw(1))); %Define k factors D = I_22*I_33 - I_23^2; ky = I_22/D; kz = I_33/D; kyz = I_23/D; % Generate panel contours Bloc = min(abs(yu-.25*c)); Floc = min(abs(yu-.75*c)); Eloc = min(abs(yl-.25*c)); Hloc = min(abs(yl-.75*c)); j=1; while abs(yu(j)-.25*c)~= Bloc j = j+1; end Bj = j; j=1; while abs(yu(j)-.75*c)~= Floc j = j+1; end Fj = j; j=1; while abs(yl(j)-.25*c)~= Eloc j = j+1; end Ej = j; j=1; while abs(yl(j)-.75*c)~= Hloc j = j+1; end Hj = j; TF_y = yu(length(yu):-1:Fj); TF_z = zu(length(zu):-1:Fj)+ .5*ts; BF_y = yu(Fj:-1:Bj); BF_z = zu(Fj:-1:Bj) + .5*ts; BE_yu = yu(Bj:-1:1); BE_yl = yl(1:Ej); BE_y = [BE_yu,BE_yl]; BE_zu = zu(Bj:-1:1) + .5*ts; BE_zl = zl(1:Ej) - .5*ts; BE_z = [BE_zu,BE_zl]; EH_y = yl(Ej:Hj); EH_z = zl(Ej:Hj) - .5*ts; HT_y = yl(Hj:length(yl)); HT_z = zl(Hj:length(zl))- .5*ts; % Generate panel lengths s_TF = pathlength(TF_y,TF_z); s_BF = pathlength(BF_y,BF_z); s_BE = pathlength(BE_y,BE_z); s_EH = pathlength(EH_y,EH_z); s_HT = pathlength(HT_y,HT_z); % Generate panel areas A_TF = ts*s_TF; A_BF = ts*s_BF; A_BE = ts*s_BE; A_EH = ts*s_EH; A_HT = ts*s_HT; % Generate panel centroids TFcent(1,1) = TF_y(1) - .5*s_TF; TFcent(2,1) = (TF_z(1) + TF_z(length(TF_z)))/2; BFcent(1,1) = BF_y(length(BF_y)) + .5*s_BF; BFcent(2,1) = BF_z(length(BF_z)) - .5*ts; BEcent(1,1) = ((BE_yu(1) + BE_yl(length(BE_yl)))/2)/2; BEcent(2,1) = ((BE_zu(1) + BE_zl(length(BE_zl)))/2); EHcent(1,1) = EH_y(1) + .5*s_EH; EHcent(2,1) = EH_z(1) + .5*ts; HTcent(1,1) = HT_y(1) + .5*s_TF; HTcent(2,1) = (HT_z(1) + HT_z(length(HT_z)))/2; % Calculate Qys and Qzs Qy_BF = A_F*(Fcent(2,1) - cw(2,1)) + A_BF*(BFcent(2,1) - cw(2,1)) + ... A_TF*(TFcent(2,1) - cw(2,1)); Qy_BE = A_F*(Fcent(2,1) - cw(2,1)) + A_BF*(BFcent(2,1) - cw(2,1)) + ... A_TF*(TFcent(2,1) - cw(2,1)) + A_B*(Bcent(2,1) - cw(2,1)) + ... A_BE*(BEcent(2,1) - cw(2,1)); Qy_EH = A_F*(Fcent(2,1) - cw(2,1)) + A_BF*(BFcent(2,1) - cw(2,1)) + ... A_TF*(TFcent(2,1) - cw(2,1)) + A_B*(Bcent(2,1) - cw(2,1)) + ... A_BE*(BEcent(2,1) - cw(2,1)) + A_E*(Ecent(2,1) - cw(2,1)) + ... A_EH*(EHcent(2,1) - cw(2,1)); Qz_BF = A_F*(Fcent(1,1) - cw(1,1)) + A_BF*(BFcent(1,1) - cw(1,1)) + ... A_TF*(TFcent(1,1) - cw(1,1)); Qz_BE = A_F*(Fcent(1,1) - cw(1,1)) + A_BF*(BFcent(1,1) - cw(1,1)) + ... A_TF*(TFcent(1,1) - cw(1,1)) + A_B*(Bcent(1,1) - cw(1,1)) + ... A_BE*(BEcent(1,1) - cw(1,1)); Qz_EH = A_F*(Fcent(1,1) - cw(1,1)) + A_BF*(BFcent(1,1) - cw(1,1)) + ... A_TF*(TFcent(1,1) - cw(1,1)) + A_B*(Bcent(1,1) - cw(1,1)) + ... A_BE*(BEcent(1,1) - cw(1,1)) + A_E*(Ecent(1,1) - cw(1,1)) + ... A_EH*(EHcent(1,1) - cw(1,1)); % Find q's fz = -(ky*V_y - kyz*V_z); fy = -(kz*V_z - kyz*V_y); qBF = fz*Qz_BF + fy*Qy_BF; qBE = fz*Qz_BE + fy*Qy_BE; qEH = fz*Qz_EH + fy*Qy_EH; qHF = 0; %due to cut at point T % Find moment separations dBF = sqrt((BFcent(2,1) - pointD(2,1))^2 + (BFcent(1,1) - pointD(1,1)^2)); dBE = sqrt((BEcent(2,1) - pointD(2,1))^2 + (BEcent(1,1) - pointD(1,1)^2)); dEH = sqrt((EHcent(2,1) - pointD(2,1))^2 + (EHcent(1,1) - pointD(1,1)^2)); % Find q_0 q_0 = (-qBF*s_BF*dBF + qBE*s_BE*dBE - qEH*s_EH*dEH)/(4*A_bar); % Generate true shear flows q_BF = qBF + q_0; q_BE = qBE + q_0; q_EH = qEH + q_0; q_HF = q_0; fprintf('\n') fprintf('The true shear flow from stringer F to B is %1.3e N/m\n',q_BF) fprintf('The true shear flow from stringer B to E is %1.3e N/m\n',q_BE) fprintf('The true shear flow from stringer E to H is %1.3e N/m\n',q_EH) fprintf('The true shear flow from stringer H to F is %1.3e N/m\n',q_HF) fprintf('\n') %***************************************** function [yu,zu,yl,zl]=airfoil(c,m,p,t,ns) % Function to generate coordinates for the upper and lower surface of an % arbitrary NACA 4-digit series airfoil. The inputs are: % % c - chord length % m - maximum camber in percentage chord % p - location of maximum camber in tenths of chord % t - maximum thickness in percentage chord % ns - desired number of contour segments % % The outputs are: % yu - y coordinates, upper contour points % zu - z coordinates, upper contour points % yl - y coordinates, lower contour points % zl - z coordinates, upper contour points % 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; %************************************** function [cent]=centroid(yu,zu,yl,zl,n) % Function to calculate the centroid of an arbitrary closed cross-section. % The inputs are: % % yu - y coordinates, upper contour % zu - z coordinates, upper contour % yl - y coordinates, lower contour % zl - z coordinates, lower contour % % The output is [cent], which is a column vector containing the y and z % coordinate of the area's centroid. asum = 0; cen_y = 0; cen_z = 0; % calculate centroid of cross-section 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,1) = cen_y/asum; cent(2,1) = cen_z/asum; %**************************************************** function sf = pathlength(y,z) sf = 0; dists = 0; L = length(y) - 1; for i = 1:L dists(i) = sqrt((y(i+1)-y(i))^2+(z(i+1)-z(i))^2); end sf = sum(dists);

Full MATLAB Code - Multi-Celled Airfoil, Shear Analysis

function vqcrew_shear_multicell(c,m,p,t,ns) close all % Define given parameters c = c+.001; ts = .002; %skin thickness tw = .003; %spar thickness tstr = .005; %stringer thickness A_B = 2e-4; %stringer areas A_E = 2e-4; A_H = 1e-4; A_F = 1e-4; A_T = A_B+A_E+A_H+A_F; %total stringer area M_y = -1250; %y direction bending moment M_z = 500; %z direction bending moment V_y = 10000; %horizontal shear force V_z = 5000; %vertical shear force %Generate airfoil contour and contour length [yu,zu,yl,zl] = airfoil(c,m,p,t,ns); su = pathlength(yu,zu); sl = pathlength(yl,zl); s = su+sl; 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; 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; 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; % Define point B stringer endpoints and stringer centroid Lbrack = .5*A_B/tstr; ay1 = min(abs(yu-.25*c-.5*tw)); by1 = min(abs(yu-.25*c+.5*tw)); ay2 = min(abs(yu-.25*c-.5*tw-.5*Lbrack)); by2 = min(abs(yu-.25*c+.5*tw+.5*Lbrack)); j=1; while abs(yu(j)-.25*c-.5*tw)~=ay1 j = j+1; end jb1 = j; j=1; while abs(yu(j)-.25*c+.5*tw)~=by1 j = j+1; end jb2 = j; j=1; while abs(yu(j)-.25*c-.5*tw-.5*Lbrack)~=ay2 j = j+1; end jb3 = j; j=1; while abs(yu(j)-.25*c+.5*tw+.5*Lbrack)~=by2 j = j+1; end jb4 = j; Bpoint1a = [yu(jb2);zu(jb2)-.5*ts]; Bpoint2a = [yu(jb4);zu(jb4)-.5*ts]; Bpoint3a = [Bpoint1a(1);Bpoint1a(2)-.5*Lbrack]; Bpoint4a = [Bpoint3a(1)-tstr;Bpoint3a(2)]; Bpoint5a = [Bpoint2a(1);Bpoint2a(2)-tstr]; Bpoint6a = [Bpoint4a(1);Bpoint5a(2)]; Bpoint1b = [yu(jb1);zu(jb1)-.5*ts]; Bpoint2b = [yu(jb3);zu(jb3)-.5*ts]; Bpoint3b = [Bpoint1b(1);Bpoint1b(2)-.5*Lbrack]; Bpoint4b = [Bpoint3b(1)+tstr;Bpoint3b(2)]; Bpoint5b = [Bpoint2b(1);Bpoint2b(2)-tstr]; Bpoint6b = [Bpoint4b(1);Bpoint5b(2)]; Bcent(1,1) = (Bpoint2a(1) + Bpoint2b(1) + Bpoint5a(1) + Bpoint5b(1) + ... Bpoint6a(1) + Bpoint6b(1) + Bpoint4a(1) + Bpoint4b(1))/8; Bcent(2,1) = (Bpoint2a(2) + Bpoint2b(2) + Bpoint5a(2) + Bpoint5b(2) + ... Bpoint6a(2) + Bpoint6b(2) + Bpoint4a(2) + Bpoint4b(2))/8; % Define point E stringer endpoints and stringer centroid Lbrack = .5*A_E/tstr; ay1 = min(abs(yl-.25*c-.5*tw)); by1 = min(abs(yl-.25*c+.5*tw)); ay2 = min(abs(yl-.25*c-.5*tw-.5*Lbrack)); by2 = min(abs(yl-.25*c+.5*tw+.5*Lbrack)); j=1; while abs(yl(j)-.25*c-.5*tw)~=ay1 j = j+1; end je1 = j; j=1; while abs(yl(j)-.25*c+.5*tw)~=by1 j = j+1; end je2 = j; j=1; while abs(yl(j)-.25*c-.5*tw-.5*Lbrack)~=ay2 j = j+1; end je3 = j; j=1; while abs(yl(j)-.25*c+.5*tw+.5*Lbrack)~=by2 j = j+1; end je4 = j; Epoint1a = [yl(je2);zl(je2)+.5*ts]; Epoint2a = [yl(je4);zl(je4)+.5*ts]; Epoint3a = [Epoint1a(1);Epoint1a(2)+.5*Lbrack]; Epoint4a = [Epoint3a(1)-tstr;Epoint3a(2)]; Epoint5a = [Epoint2a(1);Epoint2a(2)+tstr]; Epoint6a = [Epoint4a(1);Epoint5a(2)]; Epoint1b = [yl(je1);zl(je1)+.5*ts]; Epoint2b = [yl(je3);zl(je3)+.5*ts]; Epoint3b = [Epoint1b(1);Epoint1b(2)+.5*Lbrack]; Epoint4b = [Epoint3b(1)+tstr;Epoint3b(2)]; Epoint5b = [Epoint2b(1);Epoint2b(2)+tstr]; Epoint6b = [Epoint4b(1);Epoint5b(2)]; Ecent(1,1) = (Epoint2a(1) + Epoint2b(1) + Epoint5a(1) + Epoint5b(1) + ... Epoint6a(1) + Epoint6b(1) + Epoint4a(1) + Epoint4b(1))/8; Ecent(2,1) = (Epoint2a(2) + Epoint2b(2) + Epoint5a(2) + Epoint5b(2) + ... Epoint6a(2) + Epoint6b(2) + Epoint4a(2) + Epoint4b(2))/8; % Define point H stringer endpoints and stringer centroid Lbrack = .5*A_H/tstr; ay1 = min(abs(yl-.75*c-.5*tw)); by1 = min(abs(yl-.75*c+.5*tw)); ay2 = min(abs(yl-.75*c-.5*tw-.5*Lbrack)); by2 = min(abs(yl-.75*c+.5*tw+.5*Lbrack)); j=1; while abs(yl(j)-.75*c-.5*tw)~=ay1 j = j+1; end jh1 = j; j=1; while abs(yl(j)-.75*c+.5*tw)~=by1 j = j+1; end jh2 = j; j=1; while abs(yl(j)-.75*c-.5*tw-.5*Lbrack)~=ay2 j = j+1; end jh3 = j; j=1; while abs(yl(j)-.75*c+.5*tw+.5*Lbrack)~=by2 j = j+1; end jh4 = j; Hpoint1a = [yl(jh2);zl(jh2)+.5*ts]; Hpoint2a = [yl(jh4);zl(jh4)+.5*ts]; Hpoint3a = [Hpoint1a(1);Hpoint1a(2)+tstr]; Hpoint4a = [Hpoint3a(1)-tstr;Hpoint3a(2)]; Hpoint1b = [yl(jh1);zl(jh1)+.5*ts]; Hpoint2b = [yl(jh3);zl(jh3)+.5*ts]; Hpoint3b = [Hpoint1b(1);Hpoint1b(2)+tstr]; Hpoint4b = [Hpoint3b(1)+tstr;Hpoint3b(2)]; Hcent(1,1) = (Hpoint2a(1) + Hpoint2b(1) + Hpoint4a(1) + Hpoint4b(1))/4; Hcent(2,1) = (Hpoint2a(2) + Hpoint2b(2) + Hpoint4a(2) + Hpoint4b(2))/4; % Define point F stringer endpoints and stringer centroid Lbrack = .5*A_F/tstr; ay1 = min(abs(yu-.75*c-.5*tw)); by1 = min(abs(yu-.75*c+.5*tw)); ay2 = min(abs(yu-.75*c-.5*tw-.5*Lbrack)); by2 = min(abs(yu-.75*c+.5*tw+.5*Lbrack)); j=1; while abs(yu(j)-.75*c-.5*tw)~=ay1 j = j+1; end jf1 = j; j=1; while abs(yu(j)-.75*c+.5*tw)~=by1 j = j+1; end jf2 = j; j=1; while abs(yu(j)-.75*c-.5*tw-.5*Lbrack)~=ay2 j = j+1; end jf3 = j; j=1; while abs(yu(j)-.75*c+.5*tw+.5*Lbrack)~=by2 j = j+1; end jf4 = j; Fpoint1a = [yu(jf2);zu(jf2)-.5*ts]; Fpoint2a = [yu(jf4);zu(jf4)-.5*ts]; Fpoint3a = [Fpoint1a(1);Fpoint1a(2)-tstr]; Fpoint4a = [Fpoint3a(1)-tstr;Fpoint3a(2)]; Fpoint1b = [yu(jf1);zu(jf1)-.5*ts]; Fpoint2b = [yu(jf3);zu(jf3)-.5*ts]; Fpoint3b = [Fpoint1b(1);Fpoint1b(2)-tstr]; Fpoint4b = [Fpoint3b(1)+tstr;Fpoint3b(2)]; Fcent(1,1) = (Fpoint2a(1) + Fpoint2b(1) + Fpoint4a(1) + Fpoint4b(1))/4; Fcent(2,1) = (Fpoint2a(2) + Fpoint2b(2) + Fpoint4a(2) + Fpoint4b(2))/4; n = length(yu); %Define point D pointD(1,1) = (Bcent(1,1) + Ecent(1,1))/2; pointD(2,1) = (Bcent(2,1) + Ecent(2,1))/2; %Find overall stringer centroid s_cent(1,1) = (Fcent(1,1)*A_F + Bcent(1,1)*A_B + Ecent(1,1)*A_E + ... Hcent(1,1)*A_H)/A_T; s_cent(2,1) = (Fcent(2,1)*A_F + Bcent(2,1)*A_B + Ecent(2,1)*A_E + ... Hcent(2,1)*A_H)/A_T; s_cent_4str = [.2083;.0080]; % as calculated in vqcrew_4str.m % <- End Stringer Section -> % % <- Skin and Spar Section -> % % Define walls of 1/4 chord spar ylqcspar = max([Epoint1a(1),Bpoint1a(1)]); yrqcspar = min([Epoint1b(1),Bpoint1b(1)]); zlqcspar = linspace(Epoint1a(2),Bpoint1a(2),ns); zrqcspar = linspace(Epoint1b(2),Bpoint1b(2),ns); % Define walls of 3/4 chord spar yltqcspar = max([Hpoint1a(1),Fpoint1a(1)]); yrtqcspar = min([Hpoint1b(1),Fpoint1b(1)]); zltqcspar = linspace(Hpoint1a(2),Fpoint1a(2),ns); zrtqcspar = linspace(Hpoint1b(2),Fpoint1b(2),ns); % Define spar web and airfoil skin areas A_skin = s*ts; L1a = []; L1b = []; L2a = []; L2b = []; for i=1:ns-1 L1a(i) = abs(zlqcspar(i) - zlqcspar(i+1)); L1b(i) = abs(zrqcspar(i) - zrqcspar(i+1)); L2a(i) = abs(zltqcspar(i) - zltqcspar(i+1)); L2b(i) = abs(zrtqcspar(i) - zrtqcspar(i+1)); end Llqc = sum(L1a); Lrqc = sum(L1b); Lltqc = sum(L2a); Lrtqc = sum(L2b); A_qc = (Llqc+Lrqc)*.5*tw; A_tqc=(Lltqc+Lrtqc)*.5*tw; % Define airfoil skin centroid and spar web centroids qccent(1,1) = .25*c; qccent(2,1) = .5*(zlqcspar(1)+zrqcspar(1))+.25*(Llqc+Lrqc); tqccent(1,1) = .75*c; tqccent(2,1) =.5*(zltqcspar(1)+zrtqcspar(1))+.25*(Lltqc+Lrtqc); % Note: it is assumed that the airfoil skin's centroid is the same as the % general airfoil's. sk_cent = centroid(yu,zu,yl,zl,ns); % <- End Skin and Spar Section -> % % <- Calculation Section -> % % Calculate Wing Centroid A_tot = A_T + A_skin + A_qc + A_tqc; % Total wing component area cw(1,1) = (s_cent(1)*A_T + sk_cent(1)*A_skin + A_qc*qccent(1) + ... A_tqc*tqccent(1))/(A_tot); cw(2,1) = (s_cent(2)*A_T + sk_cent(2)*A_skin + A_qc*qccent(2) + ... A_tqc*tqccent(2))/(A_tot); % Calculate moment of inertia tensor components I_22 = (A_B*Bcent(2)^2 + A_B*(Bcent(1)-cw(1))^2) + ... (A_E*Ecent(2)^2 + A_E*(Ecent(1)-cw(1))^2) + ... (A_H*Hcent(2)^2 + A_H*(Hcent(1)-cw(1))^2) + ... (A_F*Fcent(2)^2 + A_F*(Fcent(1)-cw(1))^2) + ... (A_skin*sk_cent(2)^2 + A_skin*(sk_cent(1)-cw(1))^2) + ... (A_qc*qccent(2)^2 + A_qc*(qccent(1)-cw(1))^2) + ... (A_tqc*tqccent(2)^2 + A_tqc*(tqccent(1)-cw(1))^2); I_33 = (A_B*Bcent(1)^2 + A_B*(Bcent(2)-cw(2))^2) + ... (A_E*Ecent(1)^2 + A_E*(Ecent(2)-cw(2))^2) + ... (A_H*Hcent(1)^2 + A_H*(Hcent(2)-cw(2))^2) + ... (A_F*Fcent(1)^2 + A_F*(Fcent(2)-cw(2))^2) + ... (A_skin*sk_cent(1)^2 + A_skin*(sk_cent(2)-cw(2))^2) + ... (A_qc*qccent(1)^2 + A_qc*(qccent(2)-cw(2))^2) + ... (A_tqc*tqccent(1)^2 + A_tqc*(tqccent(2)-cw(2))^2); I_23 = (A_B*(Bcent(2)-cw(2))*(Bcent(1)-cw(1))) + ... (A_E*(Ecent(2)-cw(2))*(Ecent(1)-cw(1))) + ... (A_H*(Hcent(2)-cw(2))*(Hcent(1)-cw(1))) + ... (A_F*(Fcent(2)-cw(2))*(Fcent(1)-cw(1))) + ... (A_skin*(sk_cent(2)-cw(2))*(sk_cent(1)-cw(1))) + ... (A_qc*(qccent(2)-cw(2))*(qccent(1)-cw(1))) + ... (A_tqc*(tqccent(2)-cw(2))*(tqccent(1)-cw(1))); %Define k factors D = I_22*I_33 - I_23^2; ky = I_22/D; kz = I_33/D; kyz = I_23/D; % Generate panel contours Bloc = min(abs(yu-.25*c)); Floc = min(abs(yu-.75*c)); Eloc = min(abs(yl-.25*c)); Hloc = min(abs(yl-.75*c)); j=1; while abs(yu(j)-.25*c)~= Bloc j = j+1; end Bj = j; j=1; while abs(yu(j)-.75*c)~= Floc j = j+1; end Fj = j; j=1; while abs(yl(j)-.25*c)~= Eloc j = j+1; end Ej = j; j=1; while abs(yl(j)-.75*c)~= Hloc j = j+1; end Hj = j; TF_y = yu(length(yu):-1:Fj); TF_z = zu(length(zu):-1:Fj)+ .5*ts; BF_y = yu(Fj:-1:Bj); BF_z = zu(Fj:-1:Bj) + .5*ts; BE_yu = yu(Bj:-1:1); BE_yl = yl(1:Ej); BE_y = [BE_yu,BE_yl]; BE_zu = zu(Bj:-1:1) + .5*ts; BE_zl = zl(1:Ej) - .5*ts; BE_z = [BE_zu,BE_zl]; EH_y = yl(Ej:Hj); EH_z = zl(Ej:Hj) - .5*ts; HT_y = yl(Hj:length(yl)); HT_z = zl(Hj:length(zl))- .5*ts; % Generate panel lengths s_TF = pathlength(TF_y,TF_z); s_BF = pathlength(BF_y,BF_z); s_BE = pathlength(BE_y,BE_z); s_EH = pathlength(EH_y,EH_z); s_HT = pathlength(HT_y,HT_z); % Generate panel areas A_TF = ts*s_TF; A_BF = ts*s_BF; A_BE = ts*s_BE; A_EH = ts*s_EH; A_HT = ts*s_HT; % Generate panel centroids TFcent(1,1) = TF_y(1) - .5*s_TF; TFcent(2,1) = (TF_z(1) + TF_z(length(TF_z)))/2; BFcent(1,1) = BF_y(length(BF_y)) + .5*s_BF; BFcent(2,1) = BF_z(length(BF_z)) - .5*ts; BEcent(1,1) = ((BE_yu(1) + BE_yl(length(BE_yl)))/2)/2; BEcent(2,1) = ((BE_zu(1) + BE_zl(length(BE_zl)))/2); EHcent(1,1) = EH_y(1) + .5*s_EH; EHcent(2,1) = EH_z(1) + .5*ts; HTcent(1,1) = HT_y(1) + .5*s_TF; HTcent(2,1) = (HT_z(1) + HT_z(length(HT_z)))/2; % Calculate Qys and Qzs Qy_BF = A_F*(Fcent(2,1) - cw(2,1)) + A_BF*(BFcent(2,1) - cw(2,1)) + ... A_TF*(TFcent(2,1) - cw(2,1)); Qy_BE = A_F*(Fcent(2,1) - cw(2,1)) + A_BF*(BFcent(2,1) - cw(2,1)) + ... A_TF*(TFcent(2,1) - cw(2,1)) + A_B*(Bcent(2,1) - cw(2,1)) + ... A_BE*(BEcent(2,1) - cw(2,1)); Qy_EH = A_F*(Fcent(2,1) - cw(2,1)) + A_BF*(BFcent(2,1) - cw(2,1)) + ... A_TF*(TFcent(2,1) - cw(2,1)) + A_B*(Bcent(2,1) - cw(2,1)) + ... A_BE*(BEcent(2,1) - cw(2,1)) + A_E*(Ecent(2,1) - cw(2,1)) + ... A_EH*(EHcent(2,1) - cw(2,1)); Qz_BF = A_F*(Fcent(1,1) - cw(1,1)) + A_BF*(BFcent(1,1) - cw(1,1)) + ... A_TF*(TFcent(1,1) - cw(1,1)); Qz_BE = A_F*(Fcent(1,1) - cw(1,1)) + A_BF*(BFcent(1,1) - cw(1,1)) + ... A_TF*(TFcent(1,1) - cw(1,1)) + A_B*(Bcent(1,1) - cw(1,1)) + ... A_BE*(BEcent(1,1) - cw(1,1)); Qz_EH = A_F*(Fcent(1,1) - cw(1,1)) + A_BF*(BFcent(1,1) - cw(1,1)) + ... A_TF*(TFcent(1,1) - cw(1,1)) + A_B*(Bcent(1,1) - cw(1,1)) + ... A_BE*(BEcent(1,1) - cw(1,1)) + A_E*(Ecent(1,1) - cw(1,1)) + ... A_EH*(EHcent(1,1) - cw(1,1)); % Find q's fz = -(ky*V_y - kyz*V_z); fy = -(kz*V_z - kyz*V_y); qBF = fz*Qz_BF + fy*Qy_BF; qBE = fz*Qz_BE + fy*Qy_BE; qEH = fz*Qz_EH + fy*Qy_EH; qHF = 0; %due to cut at point T qqc = 0; %due to cut at point G qtqc = 0; %due to cut at point D % Find moment separations dBF = sqrt((BFcent(2,1) - pointD(2,1))^2 + (BFcent(1,1) - pointD(1,1)^2)); dBE = sqrt((BEcent(2,1) - pointD(2,1))^2 + (BEcent(1,1) - pointD(1,1)^2)); dEH = sqrt((EHcent(2,1) - pointD(2,1))^2 + (EHcent(1,1) - pointD(1,1)^2)); % Find q_0 q_0 = (-qBF*s_BF*dBF + qBE*s_BE*dBE - qEH*s_EH*dEH)/(4*A_bar); % Generate true shear flows q_BF = qBF + q_0; q_BE = qBE + q_0; q_EH = qEH + q_0; q_HF = q_0; q_qc = q_0; q_tqc = q_0; q_all = [q_BF,q_BE,q_EH,q_HF,q_qc,q_tqc]; q_all = q_all/ts; minstress = min(q_all); fprintf('\n') fprintf('The true shear flow from stringer F to B is %1.3e N/m\n',q_BF) fprintf('The true shear flow from stringer B to E is %1.3e N/m\n',q_BE) fprintf('The true shear flow from stringer E to H is %1.3e N/m\n',q_EH) fprintf('The true shear flow from stringer H to F is %1.3e N/m\n',q_HF) fprintf('The true shear flow in the c/4 spar web is %1.3e N/m\n',q_qc) fprintf('The true shear flow in the 3c/4 spar web is %1.3e N/m\n',q_tqc) fprintf('The total skin length is %1.4f m\n',s) fprintf('The minimum shear stress is %1.3f N/m^2\n',minstress) fprintf('\n') %***************************************** function [yu,zu,yl,zl]=airfoil(c,m,p,t,ns) % Function to generate coordinates for the upper and lower surface of an % arbitrary NACA 4-digit series airfoil. The inputs are: % % c - chord length % m - maximum camber in percentage chord % p - location of maximum camber in tenths of chord % t - maximum thickness in percentage chord % ns - desired number of contour segments % % The outputs are: % yu - y coordinates, upper contour points % zu - z coordinates, upper contour points % yl - y coordinates, lower contour points % zl - z coordinates, upper contour points % 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; %************************************** function [cent]=centroid(yu,zu,yl,zl,n) % Function to calculate the centroid of an arbitrary closed cross-section. % The inputs are: % % yu - y coordinates, upper contour % zu - z coordinates, upper contour % yl - y coordinates, lower contour % zl - z coordinates, lower contour % % The output is [cent], which is a column vector containing the y and z % coordinate of the area's centroid. asum = 0; cen_y = 0; cen_z = 0; % calculate centroid of cross-section 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,1) = cen_y/asum; cent(2,1) = cen_z/asum; %**************************************************** function sf = pathlength(y,z) sf = 0; dists = 0; L = length(y) - 1; for i = 1:L dists(i) = sqrt((y(i+1)-y(i))^2+(z(i+1)-z(i))^2); end sf = sum(dists);

Full MATLAB Code - Critical Buckling Shear Stresses and Buckling Modes

h = .002; v = .33; E = 72e9; b = .5; var = [.5:.1:2]; kc = ((pi^2)/32)*(9*(1+var.^2).^2)./(var.^3); D = (E*h^3)/(12*(1 - v^2)); fact = (D*pi^2)/(h*b^2); sigcr_1 = fact*kc; lambda1 = (var.^4)./(81*(1 + var.^2).^4); lambda2 = 1 + (81/625) + (81/25)*((1+var.^2)./(1+9*var.^2)).^2; lambda3 = (81/25)*((1+var.^2)./(9+var.^2)).^2; lambda = (lambda1.*(lambda2 + lambda3)).^(.5); kc_2 = (pi^2/32)*(1./(lambda.*var)); sigcr_2 = fact*kc_2; var_clamped = [1,2,inf,inf]; K_clamped = [12.7,9.5,7.38,4.31]; sig_clamped = K_clamped*(E/(1-v^2))*(h/b)^2; hold on plot(var, sigcr_1, 'b') plot(var, sigcr_2, 'r') plot(var_clamped, sig_clamped, 'g') xlabel('aspect ratio') ylabel('minimum critical buckling shear stress') hold off figure L = lambda(10); K = zeros(5,5); K(2,2) = 16*K(1,1); K(3,3) = L*(1+9*1.5^2)^2/1.5^2; K(4,4) = L*(9+1.5^2)^2/1.5^2; K(5,5) = L*(9+9*1.5^2)^2/1.5^2; K(2,1) = 4/9; K(2,3) = -4/5; K(2,4) = -4/5; K(2,5) = 36/25; K(3,2) = -4/5; K(4,2) = -4/5; K(5,2) = 36/25; Kinv = inv(K(2:5,2:5)); Cvec = Kinv*[(-4/9);0;0;0]; [x,y] = meshgrid(0:.05:1.5,0:.05:1); u_z =(sin(pi*x/.75)).*(sin(pi*y/.5)) + ... Cvec(1)*(sin(2*pi*x/.75)).*(sin(2*pi*y/.5)) + ... Cvec(2)*(sin(pi*x/.75)).*(sin(3*pi*y/.5)) + ... Cvec(3)*(sin(3*pi*x/.75)).*(sin(pi*y/.5)) + ... Cvec(4)*(sin(3*pi*x/.75)).*(sin(3*pi*y/.5)); surf(x,y,u_z) xlabel('X') ylabel('Y') zlabel('u_z')

Full MATLAB Code - Verification of Lambda

Consider the 2x2 $\mathbf {K} \;$ matrix: $\displaystyle \left[{\begin{array}{ll}{\frac {\lambda (1+\vartheta ^{2})^{2}}{\vartheta ^{2}}}&{\frac {4}{9}}\\{\frac {4}{9}}&{\frac {16\lambda (1+\vartheta ^{2})^{2}}{\vartheta ^{2}}}\end{array}}\right]$ The determinant of this matrix yields the following: $\left({\frac {\lambda (1+\vartheta ^{2})^{2}}{\vartheta ^{2}}}\right)\left({\frac {16\lambda (1+\vartheta ^{2})^{2}}{\vartheta ^{2}}}\right)-{\frac {16}{81}}=0$ Reducing, $\lambda ^{2}={\frac {\vartheta ^{4}}{9(1+\vartheta ^{2})^{2}}}$ Thus, $\lambda =\pm {\frac {1}{9}}{\frac {\vartheta ^{2}}{(1+\vartheta ^{2})}}$

Code Output - Single Celled Airfoil, Shear Analysis
>> vqcrew_shear_singlecell(c,m,p,t,ns) The true shear flow from stringer F to B is -2.065e+003 N/m The true shear flow from stringer B to E is 1.973e+003 N/m The true shear flow from stringer E to H is 1.529e+004 N/m The true shear flow from stringer H to F is 1.797e+004 N/m

Code Output - Multi-Celled Airfoil, Shear Analysis

>> vqcrew_shear_multicell(c,m,p,t,ns) The true shear flow from stringer F to B is -1.015e+003 N/m The true shear flow from stringer B to E is 3.290e+003 N/m The true shear flow from stringer E to H is 1.200e+004 N/m The true shear flow from stringer H to F is 1.374e+004 N/m The true shear flow in the c/4 spar web is 1.374e+004 N/m The true shear flow in the 3c/4 spar web is 1.374e+004 N/m The total skin length is 1.0314 m The minimum shear stress is -507299.601 N/m^2

Fig. 17 - Critical Stress vs. Aspect Ratio

Fig. 18 - Buckling Mode Shapes

Figure 17 shows the minimum critical buckling stresses versus aspect ratio. The blue line is for the 2 equation method, the red line for the 5 equation method, and the green line for clamped conditions. As expected, the clamped conditions stresses are higher. Figure 18 shows the buckling mode shapes perspective view.

The calculated minimum shear stress in the skin panels is approximately 2 orders of magnitude lower than the theoretical buckling stress. While this may be evidence of a coding error, assuming there are no coding issues this would indicate the skin panels are not particularly effective at handling shear loads.

Essays

LaTeX Code Editor

The provided LaTeX equation editor was initially quite useful, since no members of the group had any prior experience using this coding language. The interface is very straightforward and easy to use, and since the code is auto-generated, it can be cut and pasted directly into the MediaWiki editing field. The editor also included a preview feature, so it was readily apparent if any mistakes were made during the coding process. However, as the semester progressed, the experience gained from manipulating these codes in the HW reports began to render the editor less efficient. Once one has a basic working knowledge of LaTeX coding, it is far easier to just code a tilde above a bold symbol than it is to individually search out the appropriate buttons to generate the code with the editor. While it is still useful for generating larger codes, such as blank matrices and vectors, most of the required coding for the homework assignments is much more easily done by hand than it is using the editor. In short, the editor is a wonderful learning tool for LaTeX coding, but it less useful for experienced users.

InkScape

Initially, the group generated figures and charts for the homework reports using various different softwares. While there was nothing wrong with the individual submissions, compiling a report using figures generated from multiple sources is less than ideal. It is important for these sorts of assignments to maintain a consistent look and feel throughout, as you would with any sort of paper or technical report. In order to better achieve this state, the group elected to start using only InkScape to generate non-MATLAB figures. Employing InkScape allowed the group to generate consistently uniform figures for the assignments, and the difference is very noticeable in a positive way. Furthermore, InkScape is a vector graphics program, which means that a given image will appear the same regardless of of the size it is scaled to, which is very convenient for thumbnailing on MediaWiki. The end result was that the overall quality of the group's submissions was greatly improved after the adoption of InkScape as the preferred figure generating tool.

WikiEd

At the start of the semester, some of the group members elected to employ WikiEd to assist with the generation of MediaWiki content. The experience with this extension to MediaWiki closely paralleled that of the LaTeX code editor. At first, the extension was very helpful. As with the LaTeX editor, the interface for WikiEd is very straightforward and easy to use. Anytime any non-text code is required, one can simply push a button and the code is generated in the editing field. This was very helpful at the beginning of the semester when the group was still getting familiar with MediaWiki. However, as the semester went on, the extensions became more and more cumbersome. In many cases, it is much faster to simply type in the required code than it is to go hunt for the appropriate button to generate it using WikiEd. Also, it became apparent later on that if one is typing code into a WikiEd-modified editing field, and one accidentally clicks the back button on the browser (or hits the browser back key on the keyboard, if applicable), one loses all of the work entrered into the editing field since the last save. Needless to say, this is an extremely irritating bug. As with the LaTeX editor, the WikiEd extension was initially quite useful for learning purposes but became more of a burden as the semester continued.

MediaWiki vs. E-Learning

It is the general opinion of this group that E-Learning is the slightly superior option for this kind of assignment. While MediaWiki is far more powerful in terms of content handling; that is, ability to handle equations, videos, audio files, etc., it is slightly less user-friendly and has less organizational tools. Although E-Learning is notorious for server issues, generally speaking this can be avoided by approaching assignments in a timely manner. Perhaps the biggest advantage over MediaWiki when using E-Learning is that E-Learning has the built-in capability for transmitting grades to students in a confidential manner. MediaWiki has no method for handling this, and as such requires the use of email for grade distribution. All in all, while MediaWiki is a very powerful tool and a useful tool to be familiar with, E-Learning is still a superior option for class use.

Contributing Team Members

Steven Hepsworth Eas4200c.f08.vqcrew.f 21:46, 7 December 2008 (UTC)

Keith Javier Stober Eas4200c.f08.vqcrew.b 21:56, 8 December 2008 (UTC)

Kevin Klauk Eas4200c.f08.vqcrew.E 01:43, 9 December 2008 (UTC)

Casey Barnard Eas4200c.f08.vqcrew.g 04:29, 9 December 2008 (UTC)

Darin Toscano Eas4200c.f08.vqcrew.d 12:13, 9 December 2008 (UTC)

Adam Edstrand Eas4200c.f08.vqcrew.a 18:29, 9 December 2008 (UTC)

John Saxon Eas4200c.f08.vqcrew.c 18:59, 9 December 2008 (UTC)