The expression for the finite force applied to a finite unit length of a curved panel is ${\displaystyle d{\vec {F}}=q\,d{\vec {l}}}$. However, to analyze the curved panel with respects to an x and y component, we can make use of the fact that ${\displaystyle d{\vec {l}}=dl_{y}\,{\hat {j}}+dl_{z}\,{\hat {k}}}$ to render the equation:

${\displaystyle d{\vec {F}}=q\left(dl_{y}\,{\hat {j}}+dl_{z}\,{\hat {k}}\right)}$

For ${\displaystyle dl_{y}}$ and ${\displaystyle dl_{z}}$ we can implement the polar coordinate transformations:

${\displaystyle dl_{y}=dl\,cos(\theta )}$

${\displaystyle dl_{z}=dl\,sin(\theta )}$

(remembering that we are in the yz plane)

and hence:

${\displaystyle dl\,cos(\theta )=dy}$

${\displaystyle dl\,sin(\theta )=dz}$

Integrate to find the total force

${\displaystyle {\vec {F}}=\int d{\vec {F}}=q\int dy\,{\hat {j}}+dz\,{\hat {k}}}$

${\displaystyle {\vec {F}}=q\left(\Delta y{\hat {j}}+\Delta z{\hat {k}}\right)=F_{y}\,{\hat {j}}+F_{z}\,{\hat {k}}}$

For our application we can solve for the magnitude of the force finally arriving to the equation:

${\displaystyle \parallel {\vec {F}}\parallel ={\sqrt {{F_{y}}^{2}+{F_{z}}^{2}}}=q{\sqrt {{\Delta y}^{2}+{\Delta z}^{2}}}}$

We can relate our resultant force, ${\displaystyle R=\parallel {\vec {F}}\parallel }$, with the torque, ${\displaystyle T=\,2\,q\,{\bar {A}}}$ by:

${\displaystyle d{\vec {T}}={\vec {r}}\times {d{\vec {F}}}}$

by letting ${\displaystyle \rho }$ be the perpendicular distance from the point of reference to the corresponding finite force, and by remembering that ${\displaystyle d{\vec {F}}=q\,d{\vec {l}}}$, we can conclude:

${\displaystyle dT={\rho }\,q\,dl}$

${\displaystyle T=\int dT=q\int \rho \,dl=q\int 2\,dA}$

${\displaystyle T=2\,q\,{\bar {A}}}$

It is first important to define the rate of twist ${\displaystyle \theta }$ as:

${\displaystyle \theta ={\frac {d\alpha }{dx}}=rate\,\,of\,\,twist}$

AND

${\displaystyle \gamma =r{\frac {d\alpha }{dx}}=r\,\theta }$

Also, the second polar area moment of Inertia is defined as:

${\displaystyle J=\int \int _{A}r^{2}\,dA}$

To start, our equation for torque is:

${\displaystyle T=\int \int _{A}r\,\tau \,dA}$

By substituting in that ${\displaystyle \tau =G\,\gamma }$ and ${\displaystyle \gamma =r\,\theta }$ we get:

${\displaystyle T=\int \int _{A}r\,G\,(r\,\theta )=G\,\theta \int \int _{A}r^{2}dA}$

${\displaystyle T=G\,\theta \,J}$

Torsion of uniform, non-circular bars typically leads to an occurrence of warping of the cross section. Warping is when a point in the cross section under torsion, subject to deformation, will undergo a axial displacement along the length of bar. A simple example of warping is rolling up a sheet of paper then applying a torque - the rolled sheet will noticeably elongate along its axis.

Figure 4 below shows a typical depiction of displacement in a cross section view of a non-circular bar in torsion.

As can be seen above, the displacement vector, ${\displaystyle {\overrightarrow {PP'}}}$ is perpindicular to the original position vector, ${\displaystyle {\overrightarrow {OP}}}$. This method of describing the displacement vector is chosen because with the assumption that ${\displaystyle \alpha }$ is small the y and z components of the displacement vector can easily be calculated. Figure 4a below shows the y and z components of displacement for the case ${\displaystyle \beta =0}$.

By applying basic trigonometry to the figure above, ${\displaystyle u_{y}}$ and ${\displaystyle u_{z}}$ can be solved for in terms of ${\displaystyle \alpha }$:

${\displaystyle u_{y}=R(1-\cos {\alpha })\qquad u_{z}=R\sin {\alpha }}$

If the assumption that ${\displaystyle \alpha \approx 0}$ is made, then ${\displaystyle \cos {\alpha }\approx 1}$ and ${\displaystyle \sin {\alpha }\approx 0}$. These equations then become:

${\displaystyle u_{y}\approx 0\qquad u_{z}\approx {R\alpha }}$

For the general case as shown in figure 4, R becomes the initial distance OP and the y-component of the displacement of point P can be expressed as

${\displaystyle u_{y}=-(PP')sin\beta =-(OP\alpha )sin\beta =-\alpha (OP\sin \beta )=-\alpha z_{p}}$

Similarly the z-component of the displacement of point P is

${\displaystyle u_{z}=(PP')cos\beta =(OP\alpha )cos\beta =\alpha (OP\cos \beta )=\alpha y_{p}}$

Recalling the definition of ${\displaystyle \theta }$

${\displaystyle u_{y}=-\theta xz_{p}}$

${\displaystyle u_{z}=\theta xy_{p}}$

Generalizing for all points:

${\displaystyle u_{y}=-\theta xz}$ (Eq. 1)

${\displaystyle u_{z}=\theta xy}$ (Eq. 2)

For warping displacement in the axial direction (i.e. x-direction) it can be shown that

${\displaystyle u_{x}=\theta \psi (y,z)}$ (Eq. 3)

1. Kinematic Assumptions (Equations (1), (2) and (3) above) [Sun, Sec. 3.2]
2. Strain Displacement Relationship [Sun, Sec. 3.2]
3. Equilibrium Equation for Stresses [Sun, Ch. 2 , Sec. 3.2]
4. Prandtl Stress Function ${\displaystyle \phi }$ [Sun, Sec. 3.2, Eqn. 3.15]
5. Strain Compatibility Equation [Sun, Sec. 3.2, Eqn. 3.17]
6. Equation for ${\displaystyle \phi }$ [Sun, Sec. 3.2, Eqn. 3.19]
7. Boundary Conditions for ${\displaystyle \phi }$ [Sun, Sec. 3.2, Eqn. 3.24]
8. Calculation of Torque

${\displaystyle T=2\int \int _{A}^{}{\phi dA}=GJ\theta }$ [Sun, Sec. 3.2, Eqn. 3.25]

${\displaystyle J={\frac {-4}{\bigtriangledown ^{2}\phi }}\int \int _{A}^{}{\phi dA}}$

9. Formal Derivation of Torque and Shear Flow Relationship

${\displaystyle T=2q{\bar {A}}}$ [Sun, Sec. 3.5, Eqn. 3.48]

10. Calculation of Twist Angle ${\displaystyle \theta }$

${\displaystyle \theta ={\frac {1}{2G{\bar {A}}}}\oint _{}^{}{{\frac {q}{t}}ds}}$ [Sun, Sec. 3.5, Eqn. 3.56]

11. Solution of Multicell Thin-walled Sections [Sun, Sec. 3.6]

We can simplify the analysis of a complex, multi-cell airfoil by looking at a semi-circle attached to a triangle.

To find ${\displaystyle {\bar {A}}}$ we simply look at the area of each shape sepratly.

${\displaystyle {\bar {A}}={\bar {A}}_{Semi-Circle}+{\bar {A}}_{Triangle}}$

${\displaystyle {\bar {A}}={\frac {1}{2}}\left[\pi \left({\frac {b}{2}}\right)^{2}\right]+{\frac {1}{2}}\,b\,a}$

${\displaystyle {\bar {A}}=5.57}$

Remembering our derived equation for Torque, ${\displaystyle T=2\,q\,{\bar {A}}}$, we can solve for the shear flow:

${\displaystyle q={\frac {T}{2\,{\bar {A}}}}={\frac {T}{11.14}}}$

Our twist angle equation is:

${\displaystyle \theta ={\frac {1}{2\,G\,{\bar {A}}}}\oint {{\frac {q}{t}}ds}}$

However when looking at each segment of the wall, we have 3 separate sections with constant thickness, so the integral portion becomes a sum:

${\displaystyle \theta ={\frac {1}{2\,G\,{\bar {A}}}}\sum _{i=1}^{3}{\frac {q_{i}\,l_{i}}{t_{i}}}}$

${\displaystyle \theta ={\frac {q}{11.14\,G}}\left({\frac {2\,\pi }{2\,t_{1}}}+{\frac {4}{t_{2}}}+{\frac {4.47}{t_{3}}}\right)}$

${\displaystyle \theta =10{\frac {T}{G}}}$

The area of a triangle is given as:

${\displaystyle A={\frac {1}{2}}bh}$

To prove this equation, the area of triangle CDE is subtracted from the area of triangle BDE. Notice that area of a right triangle is exactly half the area of a square whose sides are equal to the legs of the triangle. Thus:

${\displaystyle A_{BDE}={\frac {1}{2}}(BD)(DE)={\frac {1}{2}}(BD)h}$

and

${\displaystyle A_{CDE}={\frac {1}{2}}(CD)(DE)={\frac {1}{2}}(CD)h}$

Subtracting the areas to solve for the triangle of interest:

${\displaystyle A_{BCE}=A_{BDE}-A_{CDE}={\frac {1}{2}}(BC-CD)h={\frac {1}{2}}bh}$

The definition of the Second Polar Area Moment of Inertia is given as:

${\displaystyle J=\int r^{2}dA}$

For a circle of radius a, the second polar moment of inertia can be solved via a double integral as

${\displaystyle J=\int r^{2}dA=\int _{0}^{2\pi }\int _{0}^{a}{r^{2}}rdrd\theta ={\frac {1}{4}}a^{4}\left[\int _{0}^{2\pi }d\theta \right]=\left[{\frac {1}{4}}a^{4}\right]\left[2\pi \right]={\frac {\pi }{2}}a^{4}}$

For a thin-walled cylinder of inner radius ${\displaystyle r_{i}=a}$ and outer radius ${\displaystyle r_{o}=b}$ the average radius can be given as:

${\displaystyle {\bar {r}}={\frac {b+a}{2}}}$

The square of the average radius is:

${\displaystyle {\bar {r}}^{2}=\left({\frac {b+a}{2}}\right)^{2}={\frac {1}{4}}\left(b^{2}+a^{2}+2ab\right)}$

If the assumption is made that the wall is thin such that thickness, ${\displaystyle t=b-a}$, is much less than a or b, then

${\displaystyle {\bar {r}}^{2}={\frac {1}{4}}\left(2b^{2}+2a^{2}\right)={\frac {1}{2}}\left(b^{2}+a^{2}\right)}$

Then,the second polar area moment of inertia for the thin-walled section can be given as

${\displaystyle J=\int r^{2}dA=\int _{0}^{2\pi }\int _{a}^{b}{r^{2}}rdrd\theta ={\frac {1}{4}}(b^{4}-a^{4})\left[\int _{0}^{2\pi }d\theta \right]=\left[{\frac {1}{4}}(b^{4}-a^{4})\right]\left[2\pi \right]={\frac {\pi }{2}}(b^{4}-a^{4})}$

Expanding the last term provides

${\displaystyle J={\frac {\pi }{2}}(b-a)(b+a)(b^{2}+a^{2})}$

Recalling the definitions for ${\displaystyle t}$, ${\displaystyle {\bar {r}}}$ and ${\displaystyle {\bar {r}}^{2}}$ from above

${\displaystyle J={\frac {\pi }{2}}t(2{\bar {r}})(2{\bar {r}}^{2})=2\pi t{\bar {r}}^{3}}$

A solid, circular cross-section has radius ${\displaystyle r_{o}^{(a)}=1cm}$. A hollow, thin-walled, circular cross-section has inner radius ${\displaystyle r_{i}^{(b)}=5cm}$ and thickness ${\displaystyle t=0.1cm}$.

a) compute ${\displaystyle A^{(a)}}$ and ${\displaystyle A^{(b)}}$

b) compute ${\displaystyle J^{(a)}}$ and ${\displaystyle J^{(b)}}$

c) compare ${\displaystyle {\frac {J^{(a)}}{J^{(b)}}}}$

d) find radius ${\displaystyle r_{i}^{(c)}}$ of thin-walled section with thickness ${\displaystyle t=0.02r_{i}^{(c)}}$ such that ${\displaystyle J^{(c)}=J^{(a)}}$ and compare ${\displaystyle {\frac {A^{(a)}}{A^{(c)}}}}$

${\displaystyle A^{(a)}=\pi (r_{o}^{(a)})^{2}=\pi (1cm^{2})=\pi cm^{2}}$

${\displaystyle A^{(b)}=\pi ((r_{i}^{(b)}+t)^{2}-(r_{i}^{(b)})^{2})=\pi ((5.1cm)^{2}-(5cm)^{2})=1.01\pi cm^{2}}$

${\displaystyle J^{(a)}={\frac {\pi }{2}}(r_{o}^{(a)})^{4}={\frac {\pi }{2}}(1cm)^{4}={\frac {\pi }{2}}cm^{4}}$

${\displaystyle J^{(b)}={\frac {\pi }{2}}((r_{i}^{(b)}t)-r_{i}^{(b)})^{4}={\frac {\pi }{2}}((5.1cm)^{4}-(5cm)^{4}=25.76\pi cm^{4}}$

${\displaystyle {\frac {J^{(a)}}{J^{(b)}}}={\frac {{\frac {\pi }{2}}cm^{4}}{25.76\pi cm^{4}}}=0.0194}$

To determine the proper ${\displaystyle r_{i}^{(c)}}$, set the second polar area moment of inertia equations for a thin-walled section and a solid section equal to one another and solve for the unkonwn:

${\displaystyle {\frac {\pi }{2}}\left[\left(1.02r_{i}^{(c)}\right)^{4}-\left(r_{i}^{(c)}\right)^{4}\right]={\frac {\pi }{2}}\left(r^{(a)}\right)^{4}={\frac {\pi }{2}}cm^{4}}$

thus,

${\displaystyle r_{i}^{(c)}=12.13cm}$

To compare areas, calculate the area ${\displaystyle A^{(c)}}$ and compare to ${\displaystyle A^{(a)}}$ determined above.

${\displaystyle A^{(c)}=\pi \left[\left(r_{o}^{(c)}\right)^{2}-\left(r_{i}^{(c)}\right)^{2}\right]=5.945\pi cm^{2}}$

Comparing the two areas,

${\displaystyle {\frac {A^{(a)}}{A^{(c)}}}=0.168}$

It is important to understand how to analyze a multi-cell section because of their widespread use in common structures such as wings. We can apply what we already know about a single cells to our analysis of multi-cells by simply adding or integrating all of the single cell sections.

If we look at the torque generated by the ${\displaystyle i^{th}}$ single cell, we get the formula:

${\displaystyle T_{i}=2\,q_{i}\,{\bar {A}}_{i}}$

To apply this equation to a multi-cell structure, we simply add all of the cells. Given that ${\displaystyle n_{c}}$ is the number of cells:

${\displaystyle T=\sum _{i=1}^{n_{c}}{T_{i}}=2\sum _{i=1}^{n_{c}}{q_{i}\,{\bar {A}}_{i}}}$

As for the rate of twist, ${\displaystyle \theta }$, since the structure is a rigid body, ${\displaystyle \theta ={\theta }_{1}=...={\theta }_{n_{c}}}$ such that the rate of twist attributed to the ${\displaystyle i^{th}}$ cell is:

${\displaystyle {\theta }_{i}={\frac {1}{2\,G_{i}\,{\bar {A}}_{i}}}\oint {{\frac {q_{i}}{t_{i}}}ds}}$

The purpose of this project is to form the basis for future structural analysis of thin walled bodies, with emphasis on airfoils. The code should be able to plot the body, determine the average area (Ā), and find the location of the centroid. These values are found through splitting the perimeter of the bodies into small discrete segments and then using properties of the cross product to determine the area from a point in space and a quadrature to determine the centroid.

Contribution from Team Aero

A code was first developed to determine the average area and centroid location of a circle and NACA 4-digit airfoil based upon user inputs. The Matlab code is posted below and Figures 8 and 9 show examples for a circle of radius 1 m and a NACA 2415 airfoil. % Matlab Script %********************************************************************* % Filename: Airfoil.m % % PURPOSE: % Calculates the average area and centroid % of a NACA airfoil or circle based on user inputs. % Plots the shape and centroid location. % % AUTHOR: % Team Aero % EAS 4200C Aerospace Structures, Fall 2008. % % Modified on: % Created on : 1 October 2008 % % DEPENDENCIES: % call: % % REMARKS: % %********************************************************************* clear % Clear system memory % Determine what shape to analyze button = questdlg('Pick which shape you wish to analyze','Shape Selection','NACA Airfoil','Circle','NACA Airfoil'); if strcmp(button, 'NACA Airfoil') % Input for variables from user prompt={ ' NACA Digit 1 (max camber in percent of chord)' ... ' NACA Digit 2 (position of max camber in tenths of chord)' ... ' NACA Digits 3-4 (maximum thicknes as percentage of chord)' ... ' Chord length in meters'... ' Y coordinate to calculate area around'... ' Z coordinate to calculate area around'... ' Number of segments to use'}; def={'2','4','15','1','0','0','40'}; % define default values of prompt parameters text_title='Please Specify Parameters'; lineNo=1; answer=inputdlg(prompt,text_title,lineNo,def); %store user input into vector m = str2num(answer{1})/100; % max camber in percent of chord p = str2num(answer{2})/10; % position of max camber in tenths of chord t = str2num(answer{3})/100; % maximum thicknes as percentage of chord cl = str2num(answer{4}); % chord length xa = str2num(answer{5}); % Y coordinate to calculate abar around ya = str2num(answer{6}); % Z coordinate to calculate abar around ns = str2num(answer{7}); % number of segments taken x1=0:1/ns:1; % mean camber line y-coords Yc(1)=0; % set initial point % Determining Points of Airfoil for i=1:round(ns*p) Yc(i+1)=m/(p^2)*(2*p*(i/ns)-(i/ns)^2); % mean camber line from y=0 to y=p theta(i+1)=atan(m/(p^2)*(2*p-2*(i/ns))); end for j=round(ns*p)+1:ns Yc(j+1)=m/((1-p)^2)*((1-2*p)+2*p*(j/ns)-(j/ns)^2); % mean camber line from y=p to y=c theta(j+1)=atan(m/((1-p)^2)*(2*p-2*(j/ns))); end for i=0:ns Yt(i+1)=t/.2*(.2969*sqrt(i/ns)-.126*(i/ns)-.3516*(i/ns)^2+.2843*(i/ns)^3-.1015*(i/ns)^4); %Airfoil thickness xu(i+1)=i/ns-Yt(i+1)*sin(theta(i+1)); yu(i+1)=Yc(i+1)+Yt(i+1)*cos(theta(i+1)); xl(i+1)=i/ns+Yt(i+1)*sin(theta(i+1)); yl(i+1)=Yc(i+1)-Yt(i+1)*cos(theta(i+1)); end % Close up trailing edge yu(ns+1)=0; yl(ns+1)=0; xu(ns+1)=1; xl(ns+1)=1; % Scale to chord line xu = cl*xu; xl = cl*xl; yu = cl*yu; yl = cl*yl; x1 = cl*x1; Yc = cl*Yc; %Average Area Calculation for i=1:ns ai=cross([xu(i)-xa yu(i)-ya 0],[xu(i+1)-xu(i) yu(i+1)-yu(i) 0]); dAu(i)=ai(3)/2; end for i=ns+1:-1:2 ai=cross([xl(i)-xa yl(i)-ya 0],[xl(i-1)-xl(i) yl(i-1)-yl(i) 0]); dAl(i)=ai(3)/2; end a_bar=abs(sum(dAu)+sum(dAl)); % Centroid Calculation for i=1:ns xbaru(i)=(xu(i+1)-xu(i))*(2*yu(i+1)+yu(i))/(3*(yu(i)+yu(i+1)))+xu(i); ybaru(i)=(yu(i)^2+yu(i+1)^2+yu(i)*yu(i+1))/(3*(yu(i)+yu(i+1))); areau(i)=(xu(i+1)-xu(i))*(yu(i)+yu(i+1))/2; xbarl(i)=(xl(i+1)-xl(i))*(2*yl(i+1)+yl(i))/(3*(yl(i)+yl(i+1)))+xl(i); ybarl(i)=(yl(i)^2+yl(i+1)^2+yl(i)*yl(i+1))/(3*(yl(i)+yl(i+1))); areal(i)=-(xl(i+1)-xl(i))*(yl(i)+yl(i+1))/2; end abar=sum(areau)+sum(areal); xareau=areau.*xbaru; yareau=areau.*ybaru; xareal=areal.*xbarl; yareal=areal.*ybarl; Centroid = [(sum(xareau)+sum(xareal))/abar, (sum(yareau)+sum(yareal))/abar]; % Plot plot(xu,yu,'k',xl,yl,'k',x1,Yc,'b',Centroid(1), Centroid(2),'+r') xlabel('Y (m)') ylabel('Z (m)') title(['NACA ' num2str(m*100) num2str(p*10) num2str(t*100) ' Airfoil']) text(.45*cl,.25*cl,['Average area = ' num2str(a_bar) ' m^2']) text(.45*cl,.2*cl,['Centroid = (' num2str(Centroid(1)) ' m, ' num2str(Centroid(2)) ' m)']) axis([0 cl -.3*cl .3*cl]); else % Input for variables from user prompt={ ' Radius in meters' ... ' Y coordinate to calculate area around'... ' Z coordinate to calculate area around'... ' Number of segments to use'}; def={'1','0','0','40'}; text_title='Please Specify Parameters'; lineNo=1; answer=inputdlg(prompt,text_title,lineNo,def); r = str2num(answer{1}); % radius xa = str2num(answer{2}); % Y coordinate to calculate abar around ya = str2num(answer{3}); % Z coordinate to calculate abar around ns = str2num(answer{4}); % number of segments taken x1=-r:2*r/ns:r; % y-coords of shape % Determines points of circle for i=1:ns+1 yu(i)=sqrt(r^2-x1(i)^2); yl(i)=-sqrt(r^2-x1(i)^2); end % Average area calculation for i=1:ns ai=cross([x1(i)-xa yu(i)-ya 0],[x1(i+1)-x1(i) yu(i+1)-yu(i) 0]); dAu(i)=ai(3)/2; end for i=ns+1:-1:2 ai=cross([x1(i)-xa yl(i)-ya 0],[x1(i-1)-x1(i) yl(i-1)-yl(i) 0]); dAl(i)=ai(3)/2; end abar=abs(sum(dAu)+sum(dAl)); % Plot plot(x1,yu,'k',x1,yl,'k',0, 0,'+r') xlabel('Y (m)') ylabel('Z (m)') title(['Circle with Radius = ' num2str(r) ' m']) text(.3*r,r,['Average area = ' num2str(abar) ' m^2']) axis([-1.1*r 1.1*r -1.1*r 1.1*r]); end Code was also made to determine the maximum number of segments needed in the calculation of the area to be within 1% accuracy of the true value for a circle and 2 successive values to be within 1% of each other for the NACA airfoil. Figures 10 and 11 shows graphs relating the effects of the number of segments on the average area. Circle code: % Matlab Script %********************************************************************* % Filename: Circle_ns_max.m % % PURPOSE: % Calculates the maximum number of segments needed to calculate the % average area of a circle in order to be within of 1% of the actual value. % Plots actual versus calculated area for different segment numbers. % % AUTHOR: % Team Aero % EAS 4200C Aerospace Structures, Fall 2008. % % Modified on: 6 October 2008 % Created on : 1 October 2008 % % DEPENDENCIES: % call: % % REMARKS: % %********************************************************************* clear % Clear memory % Initial variable values for startup r = 1; % radius ns = 1; % number of segments taken abar = 0; % area initial value j=1; % series initial value xa=0; % X coordinate to calculate abar around ya=0; % Y coordinate to calculate abar around % Iteration to find max ns while (abar < 0.99*pi*r^2) x1=-r:2*r/ns:r; % x-coords of shape for i=1:ns+1 yu(i)=sqrt(r^2-x1(i)^2); yl(i)=-sqrt(r^2-x1(i)^2); end x(j)=j; for i=1:ns ai=cross([x1(i)-xa yu(i)-ya 0],[x1(i+1)-x1(i) yu(i+1)-yu(i) 0]); dAu(i)=ai(3)/2; end for i=ns+1:-1:2 ai=cross([x1(i)-xa yl(i)-ya 0],[x1(i-1)-x1(i) yl(i-1)-yl(i) 0]); dAl(i)=ai(3)/2; end abar(j)=abs(sum(dAu)+sum(dAl)); ns=ns+1; j=j+1; end % Variables for actual area range x2=0:1:ns; area=ones(ns+1); area=area.*0.99*pi*r^2; area2=ones(ns+1); area2=area2.*1.01*pi*r^2; % Plot plot(x,abar,'k',x2,area,'b',x2,area2,'b',j-1, abar(j-1),'+r') xlabel('Number of Segments') ylabel('Computed Area (m^2)') text(10,1,'Max number of segments = 23') title('\it{Computed versus Actual Area of Circle of Radius 1 m}') NACA Airfoil Code: % Matlab Script %********************************************************************* % Filename: NACA_ns_max.m % % PURPOSE: % Calculates the maximum number of segments needed to calculate the average % area of a NACA airfoil in order to be within of 1% of the actual value. % Plots actual versus calculated area for different segment numbers. % % AUTHOR: % Team Aero % EAS 4200C Aerospace Structures, Fall 2008. % % Modified on: 6 October 2008 % Created on : 1 October 2008 % % DEPENDENCIES: % call: % % REMARKS: % %********************************************************************* clear % Clear memory m = .02; % max camber in percent of chord p = .4; % position of max camber in tenths of chord t = .15; % maximum thicknes as percentage of chord cl = 1; % chord length ns = 2; % number of segments taken j=2; % series initial value xa=0; % X coordinate to calculate abar around ya=0; % Y coordinate to calculate abar around x1=0:1/ns:1; % mean camber line x-coords Yc(1)=0; % set initial point abar(1)=0; % intial area abar(2)=.000001; % intial area % Iteration to find max ns while (abar(j)*.99 > abar(j-1)) % Determining Points of Airfoil for i=1:round(ns*p) Yc(i+1)=m/(p^2)*(2*p*(i/ns)-(i/ns)^2); % mean camber line from x=0 to x=p theta(i+1)=atan(m/(p^2)*(2*p-2*(i/ns))); end for j=round(ns*p)+1:ns Yc(j+1)=m/((1-p)^2)*((1-2*p)+2*p*(j/ns)-(j/ns)^2); % mean camber line from x=p to x=c theta(j+1)=atan(m/((1-p)^2)*(2*p-2*(j/ns))); end for i=0:ns Yt(i+1)=t/.2*(.2969*sqrt(i/ns)-.126*(i/ns)-.3516*(i/ns)^2+.2843*(i/ns)^3-.1015*(i/ns)^4); %Airfoil thickness xu(i+1)=i/ns-Yt(i+1)*sin(theta(i+1)); yu(i+1)=Yc(i+1)+Yt(i+1)*cos(theta(i+1)); xl(i+1)=i/ns+Yt(i+1)*sin(theta(i+1)); yl(i+1)=Yc(i+1)-Yt(i+1)*cos(theta(i+1)); end % Average Area Calculation j=j+1; for i=1:ns ai=cross([xu(i)-xa yu(i)-ya 0],[xu(i+1)-xu(i) yu(i+1)-yu(i) 0]); dAu(i)=ai(3)/2; end for i=ns+1:-1:2 ai=cross([xl(i)-xa yl(i)-ya 0],[xl(i-1)-xl(i) yl(i-1)-yl(i) 0]); dAl(i)=ai(3)/2; end abar(j)=abs(sum(dAu)+sum(dAl)); ns=ns+1; end x=1:j; % Coord. for x-axis % Plot plot(x,abar,'k',j, abar(j),'+r') xlabel('Number of Segments') ylabel('Computed Average Area (m^2)') text(5,.03,['Max number of segments = ' num2str(j)]) title('\it{Number of Segments Effect on Average Area of NACA 2415 Airfoil}') axis([2 j+1 0 .1]);

The following Team Aero members contributed to this report.

Jared Lee --Eas4200c.f08.aero.lee 03:06, 1 October 2008 (UTC)

Ray Strods --Eas4200c.f08.aero.strods 03:53, 7 October 2008 (UTC)

Oliver Oyama --Eas4200c.fo8.aero.oyama 17:14, 7 October 2008 (UTC)

William Kurth --Eas4200c.f08.aero.kurth 11:53, 8 October 2008 (UTC)

Gonzalo Barcia --Eas4200c.f08.aero.barcia 12:30, 8 October 2008 (UTC)