# GyroKit library and GyroRotor program

#### raghu

##### Senior Member
Bramwell does indeed move rather skillfully ( or at times what may seem carelessly) between tpp and nf axis. Obviously no one axis system is perfect and while they simplify one derivative they often complicate the other.

Bramwell's final result are all consistently in the wind axis. To arrive at the wind derivatives though he uses a combination of tpp and nf axis depending on mathematical expedience. Hence if one is trying to compare intermediate results it can get very confusing.

My sense is maybe if we compare the wind axis values ( dz/du, dx/dw... etc. )with your body axis values turned by alphaNF it may be simpler. Just a thought. No reason the current approach would not converge.

Last edited:

#### raghu

##### Senior Member
∂/∂ŵB = tan α_NF ∂/∂û - ∂/∂ŵ​

JF, in the case of hc derivatives would you not have to further adjust for the fact that tc will impact the hc derivatives as hc is in the tpp plane?

#### kolibri282

##### Super Member
Today I didn't have as much time as I had thought, the weather was fine and we made a bicycle tour to a restaurant in the country for an "al fresco" lunch but it seems we still made some progress. dcHdmu values look quite ok now as do da1dmu. Unfortunately the resolution of the 6763 charts for da1dalfa at low mu is not sufficient for a numerical result. I have backwards calculated the value of da1dalfa you would need (including sigma correction) for fd_1 and it's 0,0263. The enlarged chart of 6763 clearly shows that the value is below 0,025 so I feel your value is too high. I can not, though, claim that my value is correct. I have included only one da1dq value because here we are two orders of magnitude apart. We will start on this next, will we?
Ken Amer shows in 2309 that using a1 derivatives leads to severe errors, so we should probably use his a_prime (a') values, shouldn't we?

#### Attachments

• prog_comp_jf_jh_1-5_101115.jpg
119.8 KB · Views: 1
• alfaNf_fd1_thtaN_4pN.png
11.7 KB · Views: 1
• StabParaComp_JF_JH_all_derivatives_101015Snap1744.xls
18 KB · Views: 0
Last edited:

#### kolibri282

##### Super Member
Q1: How and where in your program do you calculate the dcTs/doM derivatives, Jean? I tried a text search but nothing turned up.
Q1 solved: found it in Rotor.java

Q2: the derivatives seem to be calculated purely numerically so I do not immediately see where the factors enter. In the analytical derivatives e.g. there is one term describing the change of cT caused by a change in mu wrt oM and one caused by a change of inflow lamda wrt to oM. Are these two taken into account in the numerical subprogram in one fell swoop and how is this done?

I had programmed the formulae of this report because they cover the case of changes in rotorspeed:
Is this what you have implemented numerically or something similar?

Last edited:

##### Newbie
My sense is maybe if we compare the wind axis values ( dz/du, dx/dw... etc. )with your body axis values turned by alphaNF it may be simpler.

Yes. One can do it directly with the GyroRotor program with body pitch set to zero. If body pitch is set to α_NF then dz/du, dx/dw... are computed in the NFP. And if body pitch is set to α_NF+a1 then dz/du, dx/dw... are computed in the TPP.

Last edited:

##### Newbie
JF, in the case of hc derivatives would you not have to further adjust for the fact that tc will impact the hc derivatives as hc is in the tpp plane?

It depends. If you want to compute the sole hc derivative, then tc doesn't play a role. If you want to compute the longitudinal moment, then you have to take into account the TPP rotation and therefore the effect of tc. That's why I prefer to work in the NFP.

By the way, there is a mistake in my equation. The good one is :

∂/∂ŵB = - tan α_NF ∂/∂û - ∂/∂ŵ​

##### Newbie
Juergen, I had an error when converting from û, ŵ to ŵB. I have updated the values. The longitudinal derivatives WRT ŵB (especially for a1) are far better.

The biggest differences are now only on high incidence values. It may be explain if you have made the assumption that rotor incidence is a small angle while I didn't.

I read the Ken Amer report 2309. What I understood is that a1 cannot be used instead of a' and it seem's logical to me because to compute a' we have to take into account the hc variation. But, it doesn't mean that a1 derivative are erroneous.

The numerical derivatives are computed in the rotor.java class with simple finite difference. The analytical derivatives are computed in the GlauertRotorDerivativeTest.java class. I did this to crosscheck derivatives. But I didn't compute the inflow derivatives analytically because it is not possible (or it is very complicated) when incidence is not a small angle.

For Ω be aware that I use normalized derivatives so that dμ/d(hat Ω) = Ω dμ/dΩ = -μ

If you look to analytical expression of tc, hc, a1 .. and if you take numerical value for dλ/dΩ, then derivatives WRT to Ω are easy to compute.

#### Attachments

• Derivatives 13 10.png
35.4 KB · Views: 1
• StabParaComp_JF_JH_all_derivatives_101315Snap1744.xls
25 KB · Views: 0
Last edited:

#### kolibri282

##### Super Member
With mu= V*cos(alfaNf)/oM dmu/doM becomes -mu/oM which is what we see in formula 10 of 2309. Your normalized value has puzzled me until a minute ago when I realized that
dμ/d(hat Ω) = Ω dμ/dΩ = -μ is probably dμ/d(hat Ω) = Ω*dμ/dΩ = -μ
in that case our oM derivatives will differ just the way the d/dmu values differ, right?

I have also calculated dcTsdoM using the formula in the report I posted in #124 and found that the two parts dcTsdmu and dcTsdlam have at least the same order of magnitude, so I don't think you can omit the latter, what do you think? ( I am using the analytical lamda derivatives also from #124). Since a lot depends on the dlamdoM values we should probably compare these first, ok?

#### Attachments

• dcTsdom_comp_101515_Snap2057.jpg
54.3 KB · Views: 0
• Formula_10_2309.PNG
6.9 KB · Views: 0
• dcTs_full_derivative.JPG
4.5 KB · Views: 0
• dcTs_full_analytical_derivative.JPG
15.9 KB · Views: 0
Last edited:

##### Newbie
in that case our oM derivatives will differ just the way the d/dmu values differ, right?
Not only ∂/∂μ but all ∂/∂Ω. My derivatives are : ∂/∂(hat Ω) = Ω * ∂/∂Ω so for instance ∂λ/∂(hat Ω) = Ω * ∂λ/∂Ω ... etc.

Here are my values of ∂λ/∂(hat Ω) :

#### Attachments

• StabParaComp_JF_JH_all_derivatives_101815.xls
26.5 KB · Views: 0
• Capture d’écran 2015-10-18 à 08.05.38.png
4.2 KB · Views: 0

#### kolibri282

##### Super Member
My dlamdoM values show a similar trend as yours but are consistently smaller and even more so at lower mu. I have calculated dcTsdoM using only the -mu*dCtsdmu term (#1) and using the two parts of the partial dcTsdoM derivative (#2, the lower formula in #128). Since lamda decreases with increasing mu the difference becomes smaller but does not seem negligible to me.

#### Attachments

• dcTsdom_comp_101815_Snap1957.jpg
21.1 KB · Views: 0
• StabParaComp_JF_JH_all_derivatives_101815Snap2019.xls
19.5 KB · Views: 0
Last edited:

##### Newbie
Juergen, you cannot write in our context :

∂Ct/∂Ω = ∂Ct/∂μ * ∂μ/∂Ω + ∂Ct/∂λ * ∂λ/∂Ω​

because λ is not an independent variable and depends on μ. When computing partial derivative you must take care of your independent variables and stay with them.

It is correct to write : Ct=Ct(μ,λ)=Ct(V,α,Ω) and then ∂Ct/∂Ω = ∂Ct/∂μ * ∂μ/∂Ω + ∂Ct/∂λ * ∂λ/∂Ω, but in this context, the value ∂Ct/∂μ is computed with λ fixed. This is not we have done because our value of ∂Ct/∂μ take into account the ∂λ/∂μ.

It is more clear if we take the analytic expression of Ct : Ct=a/4 [ 2/3 (θ0 + 3/2μ^2) + λ].

Our value of ∂Ct/∂μ is ∂Ct/∂μ = a/4 (2 μ θ0 + ∂λ/∂μ). But if you want to use your first formula you should have write : ∂Ct/∂μ = a/4 (2 μ θ0). With this value you can write :

∂Ct/∂Ω = ∂Ct/∂μ * ∂μ/∂Ω + ∂Ct/∂λ * ∂λ/∂Ω
∂Ct/∂Ω = a/4 (2 μ θ0) * -μ + a/4 * ∂λ/∂Ω
∂Ct/∂Ω = a/4 [-2 μ^2 θ0 + ∂λ/∂Ω]​

The value I put in the Excel file includes ∂λ/∂Ω.

Last edited:

#### kolibri282

##### Super Member
Over the last weeks I didn't have much time and so didn't make much progress. I have programmed another formula for dcHdmu which is from the Princeton University Department of Aerospace Report (DRA) 695. It is very close to your value for fd1 but deviates markedly for higher mu values. The second froumla is from Bramwell which is also used only slightly modified by Gessow and Meyers. Here the values for higher mu are closer but for low mu are far apart. I will leave the dcHdmu values for now because I have no idea how to further verify them. One small bit of information is the graph of dcHdmu from naca-2655. The measured values (solid line) drop off very quickly for low mu and then level out quickly towards higher mu values. This seems to indicate that the DRA-695 values my be better but of cause this is no proof.

Thanks for the thorough explanation regarding the derivatives with respect to omega. After pondering it for quite a while I feel I have still not understood it completely and have two questions:

a) I think that ∂μ/∂Ω is equal to -μ/Ω then:
∂Ct/∂Ω = a/4 (2 μ θ0) * (-μ/Ω) + a/4 * ∂λ/∂Ω
∂Ct/∂Ω = a/4 [-2 μ^2/Ω * θ0 + ∂λ/∂Ω]

b) You state that your values include ∂λ/∂Ω. Is that because of the numerical procedure you use? If not I actually don't see where the relationship between μ and λ enters in the process.

#### Attachments

• dcHdmu_comp_1108.JPG
12.1 KB · Views: 0
• dcHdmu_comp_2655_1108.JPG
30.6 KB · Views: 0
Last edited:

##### Newbie
I should have write "hat Ω" because ∂μ/∂(hat Ω) = -μ. I suggest that Ω will be understood as "hat Ω" in future post (I don't know the ascii code to write Ω with a hat like ŵ for example).

Yes ∂λ/∂Ω is numerically computed. If you look to the the java class "GlauertRotorDerivativeTest.java", lines 228 to 232 are the numerical computation of λ derivatives.

I suggest that we to try to compare ∂Q/∂Ω.

#### kolibri282

##### Super Member
So far I have found two reports containing ∂Q/∂Ω, below are the two formulae. The first one, from naca 2655, does contain a term ∂v/∂V (where lower case v is the induced velocity in the rotor disk), which is another way of writing ∂λi/∂μ, adopting your (hat Ω) then -μ*∂λi/∂μ would be ∂λi/∂Ω and to me it looks as if the whole formula is a total derivative with respect to Ω, so it should give values we could compare to your values, right? Please let me know if I am wrong since then I would wast my time. Otherwise I will calculate the values for our five example μ values and will hopefully be able to present them tomorrow evening. I will also program the second formula which is from naca 2309. The formula states that it is for the case of constant disk angle of attack. Since this value is not varied in the formula of naca-2655 the two formulae should give the same results, or am I missing something here? Generally ∂Q/∂Ω should be ∂Q/∂μ*∂μ/dΩ assuming that ∂Q/∂μ is the total derivative with respect to μ, right?

Thanks for your support and have a nice weekend!

Cheers,

Juergen

#### Attachments

• dcQdmu_2309.png
16.2 KB · Views: 1
• dcQdmu_2655.png
27.7 KB · Views: 1
Last edited:

#### kolibri282

##### Super Member
For about two weeks I have tried now to make some headway with the dcQ/doM values. If we use Jean's oM_hat these are simply the dcQ/dmu values scaled by -mu, so I tried to verfy dcQ/dmu which succeeded fairly well, apart from one annoying difference (see below). First I switched to helicopter mode in my program and calculated values for the UH-1 model at mu=0.3. A dcQ/dmu value of about -0.24 was obtained with the 2655 and the 2309 formula. This value corresponded well to the 6763 diagram for mu=0.3. Unfortunately the 6763 forumla was off by almost an order of magnitude. Close inspection revealed a second error in the formula, delta2 is missing with the dt56/dmu value, but this didn't help much. The 6763 formula is actually the derivative of the decellerating (equation 11) and accellerating torques (equation 9) from naca-716. It is straightforward to derive this formula with respect to mu (that's how I found the missing delta2 error) but the values obtained are not in accordance with diagram 7.5-24 (p. 219) of 6763. My next idea was, that rotor solidity sigma might have a larger influence and so I used my 67-63 helicopter model, which has a solidity of 0.1 just like it is used in the report. At mu=0.3 the 2655 and 6763 values were close to within about 15% but they were only about half the value of diagram 7.5-24. I will skip the rest of the nitty gritty of what I did and found because finally I think we can at least say that the dcQdmu values are in the ball park and since my values using the formula from naca-2655 match your values quite well, Jean, I would like to leave it at that for the moment. The dcQ/doM_hat values are below and in the attached spreadsheet. I think we should move on, what do we tackle next?

#### Attachments

• dcQdoM_hat.JPG
24.7 KB · Views: 0
• dcQdoM_6763.png
22 KB · Views: 0
• naca-2309_allDerivatives.png
20.1 KB · Views: 0
• StabParaComp_JF_JH_all_derivatives_including_dcQdoM_hat_112215Snap1808.xls
20.5 KB · Views: 0
Last edited:

##### Newbie
Juergen , sorry for the delay in my response. I have not had much time to work on this in the past days. My expression of dcQ/doM is complex. You can find it in the document I just finish and which gives analytic expression of all my derivatives. You can download it here.

I suggest now that we move on the computation of the stability derivative matrix which allows to compute directly the stability mode. I will give you by the end of the week the values I got.

##### Newbie
Juergen, here is a case of stability computation :

Configuration data :

Rotor : b(number of blade)=2, R(radius)=4.2m, c(chord)=0.2m, Ib(blade flapping m.o.i.)=100 kg.m^2, Ir(rotor lag m.o.i)=200 kg.m^2, θTW(blade twist)=0deg, a(lift curve slope)=5.7, δ(profil drag)=0.013, B(tip loss factor)=0.976, xS=-0.16m, zS=-1.5m, θ0(blade pitch)=4.5deg

Fuselage : SF(front aera)=1m^2, CDF(drag coefficient)=1, aF(lift curve slope)=0, xF=0, zF=0

Tail plane : SH(aera)=0.9m^2, CDH(drag coefficient)=0, aH(lift curve slope)=3.5, ΦH(tilt)=0deg, xH=-1.6m, zH=0.5m

Propeller : Rp(radius)=0.886m, ΦP(tilt)=0deg, xP=-0.9m, zP=-0.2m

Gyro : Iyy(longitudinal m.o.i)=500 kg.m^2, m(mass)=400kg

x, y are location in body-coordinates system relative to CG. In this example I didn't take effect of rotor wake and propeller wake at tail plane.

Trim :

At 32 m/s, I get : Ω (rotor speed)=341.05rpm, αF(fuselage incidence)=-0.54deg, αR(rotor incidence)=3.91deg, TP(propeller thrust)=114.34daN, Fuselage drag=62.72daN, Rotor thrust=397.95daN, a1(longitudinal flapping)=2.90deg, tail lift=-1.86daN.

Derivatives :

Rotor derivatives :
Xu -0,043 Xw -0,162 Xq 1,365 Xθ 0,000 XΩ -0,041
Zu -0,135 Zw -0,961 Zq -0,068 Zθ 0,000 ZΩ -0,440
Mu 0,034 Mw 0,072 Mq -1,647 Mθ 0,000 MΩ -0,007
Qu 0,138 Qw 0,719 Qq -4,240 Qθ 0,000 QΩ -0,118

Tail plane derivatives :
Xu -0,000 Xw -0,003 Xq -0,005 Xθ 0,000 XΩ 0,000
Zu 0,001 Zw -0,154 Zq -0,246 Zθ 0,000 ZΩ 0,000
Mu 0,002 Mw -0,199 Mq -0,317 Mθ 0,000 MΩ 0,000
Qu 0,000 Qw 0,000 Qq 0,000 Qθ 0,000 QΩ 0,000

Fuselage derivatives :
Xu -0,098 Xw 0,000 Xq 0,000 Xθ 0,000 XΩ 0,000
Zu 0,000 Zw -0,049 Zq 0,000 Zθ 0,000 ZΩ 0,000
Mu -0,000 Mw 0,000 Mq 0,000 Mθ 0,000 MΩ 0,000
Qu 0,000 Qw 0,000 Qq 0,000 Qθ 0,000 QΩ 0,000

Total derivatives :
Xu -0,141 Xw -0,165 Xq 1,361 Xθ -9,810 XΩ -0,041
Zu -0,133 Zw -1,164 Zq 31,686 Zθ 0,093 ZΩ -0,440
Mu 0,036 Mw -0,127 Mq -1,964 Mθ 0,000 MΩ -0,007
Qu 0,138 Qw 0,719 Qq -4,240 Qθ 0,000 QΩ -0,118

Modes :

Short period : Period=3,25s Time to half=0,47s Roots : -1,462 (1/s) ±i 0,308 (cycles/s)

Long period : Period=19,50s Time to half=24,97s Roots : -0,028(1/s) ±i 0,051(cycles/s)

Rotor aperiodic : Time to half=1,71 s Roots : -0,406 (1/sec)

All modes are stables. Ask questions if something is unclear.

#### Jean Claude

##### Junior Member
Jean, It seems to me that your moments are not balanced.
Rotor : 133 mN nose up
Tail : 30 mN nose up
Prop. : 228 mN nose down

##### Newbie
Jean-Claude,

our differences come from the rotor rear force H. I got : rotor thrust T=3979.6N and rear force H=244.9N. This leads to an angle of 3.52deg of the total rotor force relative to the no-feathering plane (the total rotor force is not perfectly perpendicular to the the tip-path plane).

The angle to take into account is then 3.91+3.52+0.54=7.97deg instead of 7.35deg. The corresponding moment is 198 m.N.

#### Jean Claude

##### Junior Member
You certainly have right, Jean. Thank you for that clarification.

Replies
11
Views
649
Replies
12
Views
469
Replies
6
Views
585
Replies
2
Views
320
Replies
13
Views
523