Interactive “Snow Cone” airfoil shear flow

by. Aaron S. Holben

Calculating the centroid!

For ease of reference we’re going to calculate the centroid from the far left. First lets calculate the x centroid. Its obvious where the y centroid’s distance will be so we will skip that and pre assign it. It is not mentioned in the image but the length of the smaller spar is R.

Before going for the centroid calculations we need to get the actual lengths of every cell piece for later.

Length Calculations

BigHypotenuse = (Sqrt[(R^2)+(((4R))^2)])

Angle of each side of the big triangle.

α = ArcTan[(R/(4R))]//N

Angle  = ArcSin[((4R)/BigHypotenuse)]((180/π))//N

Cell3Hypotenuse = (Sqrt[(((((1/2))R))^2)+ (((2R))^2)])

Ls1 = π R

Ls2 = 2R

Ls3 = (Sqrt[17])R - .5(Sqrt[17])R

Ls4 = .5(Sqrt[17])R

Ls5 = R

Ls6 = .5(Sqrt[17])R

Ls7 = .5(Sqrt[17])R

Now that has been completed we can go ahead and calculate the centroid of the system

Centroid Calculation

CentroidX = ((((R-((2R)/π))(π R t))+R(2 R t) + 2(3R)((Sqrt[17])R t)+3R(R t))/(π R t +2R t + 2(Sqrt[17])R t+R t ))//N

Due to Symmetry

CentroidY = R

Calculate the Bending Stiffnesses

Now we will calculate (Subscript[H, 22]) and (Subscript[H, 33]) but there will be no need to calculate (Subscript[H, 23]) since the section is symmetrical along one axis

(Subscript[H, 23])=0;

H22= Ε((∫_0^π)(((-R Cos[θ]))^2)t R\[DifferentialD]θ+(1/12)t(((2R))^3)+(1/12)t(((R))^3)+2((1/12)((((Sqrt[17])R))^3)t ((Sin[α])^2))+(Sqrt[17])R t((((R/2)))^2))

This calculation will possibly be used for additional flow in the X direction.

c1 = CentroidX - R

c2 = 3R -CentroidX

Bear in mind that a Stiffness calculation that invokes using (t^3)will be counted as practically zero due to the fact that the moment of innertia contribution of t on any other higher power is next to nothing. So both spar base calculations are 0.

H33 = Ε((∫_0^π)(((-c1 - R Sin[θ]))^2)t R\[DifferentialD]θ+ (0+ (2R t)(((-c1))^2))+ (0+ (R t)(((c2))^2))+2((1/12)((((Sqrt[17])R))^3)t ((Cos[α])^2)+((Sqrt[17])R t(((2R))^2))))//N

Find Open Shear flow for each section

Now we will use an integral to find Q the stiffniss static moment for every section of the airfoil

Q21 = Ε(∫_0^θ)(-R Cos[θ])R t\[DifferentialD]θ

Q31 = Ε(∫_0^θ)(-c1 + R Sin[θ])R t\[DifferentialD]θ

Q22 = Ε(∫_0^s2)(-R+s2)t\[DifferentialD]s2//Simplify

Q32 = Ε(∫_0^s2)(-c1)t\[DifferentialD]s2

Q23 = Ε(∫_0^s3)(-R+s3 Sin[α])t\[DifferentialD]s3//Simplify

Q33 = Ε(∫_0^s3)(-c1+s3 Cos[α])t\[DifferentialD]s3

Q24 =Ε(∫_0^s4)(-((R/2))+s4 Sin[α])t\[DifferentialD]s4//Simplify

Q34 = Ε(∫_0^s4)(c2+s3 Cos[α])t\[DifferentialD]s4

Q25 = Ε(∫_0^s5)(((-R)/2)+s5)t\[DifferentialD]s5

Q35 = Ε(∫_0^s5)(c2)t\[DifferentialD]s5

Q26 = Ε(∫_0^s6)((R/2)+s6 Sin[α])t\[DifferentialD]s6

Q36 = Ε(∫_0^s6)(c2+s6 Cos[α])t\[DifferentialD]s6

Q27 = Ε(∫_0^s7)(s7 Sin[α])t\[DifferentialD]s7

Q37 = Ε(∫_0^s7)((5R-CentroidX)+s7 Cos[α])t\[DifferentialD]s7

Calculate “Open” Shear flow for the airfoil system

The full equation for each flow calculation is listed as:

f(s) = c+(((Subscript[Q, 3])(s)(Subscript[H, 23])-(Subscript[Q, 2])(s)(Subscript[H, 33]))/(ΔH))(Subscript[V, 3])-(((Subscript[Q, 3])(s)(Subscript[H, 22])-(Subscript[Q, 2])(s)(Subscript[H, 23]))/(ΔH))(Subscript[V, 2])

Where... ΔH = (Subscript[H, 22])(Subscript[H, 33])

We can reduce the equation a little bit since our structure is symmetrical implying our (Subscript[H, 23]) is 0. The equation is now...

f(s)= -(((Subscript[Q, 2])(s))/(Subscript[H, 22]))(Subscript[V, 3])-(((Subscript[Q, 3])(s))/(Subscript[H, 33]))(Subscript[V, 2])

Now we can calculate the open flow equations!

fo1 =(Q21/H22)V3 - (Q31/H33)V2

fo2 =  -(Q22/H22)V3 - (Q32/H33)V2

fo3 = -(Q23/H22)V3 - (Q33/H33)V2

fo4 = fo3 - (Q24/H22)V3-(Q34/H33)V2/.s3->Ls3

fo5 = fo3 -(Q25/H22)V3-(Q35/H33)V2/.s3->Ls3

This required a second step since theres Flow s5 and s7 going into s6.

fo6p = fo7 -(Q26/H22)V3-(Q36/H33)V2/.s7->Ls7

fo6 = fo6p + fo5/.s5->Ls5

fo7 = fo4 - (Q27/H22)V3-(Q37/H33)V2/.s4->Ls4

Calculate the “closed” additional flow for each flow equation

In this section we will calculate the closed additional flow that makes the flow linear rather than first order. (Which means its not a cut tube basically).

warp1 = (∫_0^(2R))((fo2 +fc1-fc2)/(G t))\[DifferentialD]s2-(∫_0^π)((fo1 -fc1)/(G t))R\[DifferentialD]θ //Simplify

warp2 = (∫_0^Ls3)((fo3+fc2)/(G t))\[DifferentialD]s3+(∫_0^Ls5)((fo5+fc2-fc3)/(G t))\[DifferentialD]s5+(∫_0^Ls6)((fo6+fc2)/(G t))\[DifferentialD]s6-(∫_0^Ls2)((fo2+fc1 -fc2)/(G t))\[DifferentialD]s2//Simplify

warp3 = (∫_0^Ls4)((fo4+fc3)/(G t))\[DifferentialD]s4+(∫_0^Ls7)((fo7+fc3)/(G t ))\[DifferentialD]s7-(∫_0^Ls5)((fo5 +fc2-fc3)/(G t))\[DifferentialD]s5

Clear[Soln1]

Soln1 = Solve[{warp1==0,warp2==0,warp3 ==0},{fc1,fc2,fc3}][[1]]//FullSimplify

{fc1,fc2,fc3}={fc1,fc2,fc3}/.Soln1

Calculate the Net shear flow equations (without added V loads)

Now all we need to do is add the flows to the respective open flow equations!

f1 = fo1 - fc1

f2 = fo2 +fc1

f3 = fo3+fc2

f4 = fo4+fc3

f5 = fo5+fc2-fc3

f6 = fo6 +fc2

f7 = fo7+fc3

Find the shear center of the system!

For this calculation we will find where the shear center is of the system. The shear center is where the forces V3 and V2 need to be so that NO ADDITIONAL SHEAR FLOW IS ADDED TO THE SYSTEM. We need to know this location to do our interactive model. It is heavily implied that the V2 shear center is along the axis of symmetry on i2. The only shear center we need to calculate is the V3 one.

Mcalc = -(∫_0^π)R f1 R\[DifferentialD]θ+R(∫_0^Ls3)f3 Cos[α]\[DifferentialD]s3+((R/2))(∫_0^Ls4)f4 Cos[α]\[DifferentialD]s4 +(c2)(∫_0^Ls4)f4 Sin[α]\[DifferentialD]s4+(c2)(∫_0^Ls5)f5\[DifferentialD]s5+((R/2))(∫_0^Ls6)f6 Cos[α]\[DifferentialD]s6+(c2)(∫_0^Ls6)f6 Sin[α]\[DifferentialD]s6+(4R)(∫_0^Ls7)f7 Sin[α]\[DifferentialD]s7

Mc2 = Mcalc/.V2->0

Solv1 = Solve[Mc2+d V3==0,d][[1]]//Simplify

Clear[d]

{d}={d}/.Solv1

Where d is the distance from the first spar in the airfoil.

Set up the interactive model!

The additional flow is simple to add on into the equation now that we know the shear center. It is simple as:

fadditional = (Moment/(2Area of Cell))

All we need to do is appropiately calculate a moment for the whole system accordingly. I will allow the shear force to only reach the ends of the whole section and not off of it. That was the purpose of calculating the shear center. Realistically the forces can only directly act on the section right?

AC1 = ((π (R^2))/2);

AC2 = .5(2R)4R -.5(R)2R;

AC3 = .5(R)2R;

Set up the twist equations along with the additional moment equations to “add” into the shear stress graph formulas

eqn1 = T == 2AC1 ft1 +2 AC2 ft2 +2AC3 ft3;

eqn2 = κ1 == (1/(2AC1))(((ft1 π R)/(G t))+(((ft1 -ft2)2R)/(G t)))//Simplify

eqn3 = κ1 ==(1/(2AC2))(((ft2 Ls3)/(G t))+(((ft2-ft3)Ls5)/(G t))+((ft2 Ls6)/(G t))+(((-ft1+ft2)2R)/(G t)))//Simplify

eqn4 = κ1 ==(1/(2 AC3))(((ft3 Ls4)/(G t))+((ft3 Ls7)/(G t))+(((ft3 -ft2)Ls5)/(G t)))//Simplify

Solv2 = Solve[{eqn1,eqn2,eqn3,eqn4},{ft1,ft2,ft3,κ1}]

Clear[ft1]

Clear[ft2]

Clear[ft3]

ft1 = ft1/.Solv2[[1]]//Simplify

ft2 = ((0.10210358028582928` T)/(R^2));

ft3 = ((0.052132251746094774` T)/(R^2));

Clear[T]

We’re going to assign V2 and V3 to a common force multiplier of P for the sake of the graphing simulation. We can adjust the magnitude of these forces in relation to each other later for analysis use.

V2 = v2 P;
V3 = v3 P;

T = V3(x R) +V2(y R);

fa1 = f1 -ft1//FullSimplify

fa2 = f2 +ft1-ft2//FullSimplify

fa3 = f3 + ft2

fa4 = f4 + ft3

fa5 = f5 +ft2-ft3

fa6 = f6 + ft2

fa7 = f7 +ft3

Clear[F1,F2,F3,F4,F5,F6,F7]

Change the graphs to be unitless in terms of force P

{F1,F2,F3,F4,F5,F6,F7}={((fa1 R)/P),((fa2 R)/P),((fa3 R)/P),((fa4 R)/P),((fa5 R)/P),((fa6 R)/P),((fa7 R)/P)}/.{s2->η Ls2,s3->η Ls3, s4->η Ls4, s5->η Ls5, s6->η Ls6, s7->η Ls7}//Simplify

{F1s,F2s,F3s,F4s,F5s,F6s,F7s}={((f1 R)/P),((f2 R)/P),((f3 R)/P),((f4 R)/P),((f5 R)/P),((f6 R)/P),((f7 R)/P)}/.{s2->η Ls2,s3->η Ls3, s4->η Ls4, s5->η Ls5, s6->η Ls6, s7->η Ls7}//Simplify

This section is just a regurgitation of what was stated earlier for our graphs. FXs graphs represent the original trend (no additional twist or torque added by shear forces) and FX is graphs that have the additional shear flow element. Notice on the interactive model when the x and y position is set to their shear centers we get the same amount of shear stress.

F1s

F2s

F3s

F4s

F5s

F6s

F7s

F1

F2

F3

F4

F5

F6

F7

Leftcalc = R+d

Rightcalc = 4R-d

Manipulate[Row[{
Plot[{-0.06034947630082764` v2+0.2839618757131349` v3+0.019893628989802483` v2 θ+0.017354978191817746` v2 Cos[θ]-0.24760093240965586` v3 Sin[θ],v3 (0.2839618757131349` -0.09011798982574327` x)+v2 (-0.06034947630082764`-0.09011798982574327` y+0.019893628989802483` θ)+0.017354978191817746` v2 Cos[θ]-0.24760093240965586` v3 Sin[θ]},{θ,0,π},ImageSize->300,Frame->True,FrameLabel->{"Θ","Shear","Nose",""}],
Plot[{v2 (0.04299449810900989` +0.039787257979604966` η)+v3 (-0.2839618757131349`+0.4952018648193117` η-0.4952018648193117` (η^2)),v2 (0.04299449810900989` -0.011985590460086001` y+0.039787257979604966` η)+v3 (-0.2839618757131349`-0.011985590460086001` x+0.4952018648193117` η-0.4952018648193117` (η^2))},{η,0,1},ImageSize->300,Frame->True,FrameLabel->{"η","Shear","First Spar",""}],
Plot[{v3 (-0.3998735704869084`+0.5104423986632151` η-0.12761059966580376` (η^2))+v2 (0.10859930073829971` +0.04101176680090259` η-0.03577820410757778` (η^2)),v3 (-0.3998735704869084`+0.10210358028582928` x+0.5104423986632151` η-0.12761059966580376` (η^2))+v2 (0.10859930073829971` +0.10210358028582928` y+0.04101176680090259` η-0.03577820410757778` (η^2))},{η,0,1},ImageSize->300,Frame->True,FrameLabel->{"η","Shear","Bottom Middle Section",""}],
Plot[{v2 (0.10876564614542962` -0.10210104962940854` η)+v3 (0.006959733717652572` +0.25522119933160753` η-0.12761059966580376` (η^2)),v2 (0.10876564614542962` +0.052132251746094774` y-0.10210104962940854` η)+v3 (0.006959733717652572` +0.052132251746094774` x+0.25522119933160753` η-0.12761059966580376` (η^2))},{η,0,1},ImageSize->300,Frame->True,FrameLabel->{"η","Shear","Bottom End Section",""}],
Plot[{v2 (0.010300779979519711` -0.014816327393833009` η)+v3 (0.35883029379026166` +0.12380046620482793` η-0.12380046620482793` (η^2)),v2 (0.010300779979519715` +0.0499713285397345` y-0.014816327393833009` η)+v3 (0.35883029379026166` +0.0499713285397345` x+0.12380046620482793` η-0.12380046620482793` (η^2))},{η,0,1},
ImageSize->300,Frame->True,FrameLabel->{"η","Shear","Second Spar",""}],
Plot[{v3 (0.3657900275079141` -0.25522119933160753` η-0.12761059966580376` (η^2))+v2 (-0.13573020463527855`-0.030544641414252963` η-0.03577820410757778` (η^2)),v3 (0.3657900275079141` +0.10210358028582928` x-0.25522119933160753` η-0.12761059966580376` (η^2))+v2 (-0.13573020463527855`+0.10210358028582928` y-0.030544641414252963` η-0.03577820410757778` (η^2))},{η,0,1},
ImageSize->300,Frame->True,FrameLabel->{"η","Shear","Top Middle Section",""}],
Plot[{v3 (0.1345703333834563` -0.12761059966580376` (η^2))+v2 (0.006664596516021082` -0.10210104962940852` η-0.03577820410757778` (η^2)),v3 (0.1345703333834563` +0.052132251746094774` x-0.12761059966580376` (η^2))+v2 (0.006664596516021082` +0.052132251746094774` y-0.10210104962940852` η-0.03577820410757778` (η^2))},{η,0,1},ImageSize->300,Frame->True,FrameLabel->{"η","Shear","Top End Section",""}]}],{x,-0.8527719802570544,4.147228019742945},{y,-1,1},{v2,1,10},{v3,1,10}]