1 // SIS Twin Otter Iced aircraft Nonlinear model
2 // Version 020409
3 // read readme_020212.doc for information
4
5 // 11-21-2002 (RD) Added e code from Kishwar to fix negative lift problem at
6 // high etas
7
8 #include "uiuc_iced_nonlin.h"
9
Calc_Iced_Forces()10 void Calc_Iced_Forces()
11 {
12 // alpha in deg
13 double alpha;
14 double de;
15 double eta_ref_wing = 0.08; // eta of iced data used for curve fit
16 double eta_ref_tail = 0.20; //changed from 0.12 10-23-2002
17 double eta_wing;
18 double e;
19 //double delta_CL; // CL_clean - CL_iced;
20 //double delta_CD; // CD_clean - CD_iced;
21 //double delta_Cm; // CM_clean - CM_iced;
22 double delta_Cm_a; // (Cm_clean - Cm_iced) as a function of AoA;
23 double delta_Cm_de; // (Cm_clean - Cm_iced) as a function of de;
24 double delta_Ch_a;
25 double delta_Ch_e;
26 double KCL;
27 double KCD;
28 double KCm_alpha;
29 double KCm_de;
30 double KCh;
31 double CL_diff;
32 double CD_diff;
33
34
35
36 alpha = Std_Alpha*RAD_TO_DEG;
37 de = elevator*RAD_TO_DEG;
38 // lift fits
39 if (alpha < 16)
40 {
41 delta_CL = (0.088449 + 0.004836*alpha - 0.0005459*alpha*alpha +
42 4.0859e-5*pow(alpha,3));
43 }
44 else
45 {
46 delta_CL = (-11.838 + 1.6861*alpha - 0.076707*alpha*alpha +
47 0.001142*pow(alpha,3));
48 }
49 KCL = -delta_CL/eta_ref_wing;
50 eta_wing = 0.5*(eta_wing_left + eta_wing_right);
51 if (eta_wing <= 0.1)
52 {
53 e = eta_wing;
54 }
55 else
56 {
57 e = -0.3297*exp(-5*eta_wing)+0.3;
58 }
59 delta_CL = e*KCL;
60
61
62 // drag fit
63 delta_CD = (-0.0089 + 0.001578*alpha - 0.00046253*pow(alpha,2) +
64 -4.7511e-5*pow(alpha,3) + 2.3384e-6*pow(alpha,4));
65 KCD = -delta_CD/eta_ref_wing;
66 delta_CD = eta_wing*KCD;
67
68 // pitching moment fit
69 delta_Cm_a = (-0.01892 - 0.0056476*alpha + 1.0205e-5*pow(alpha,2)
70 - 4.0692e-5*pow(alpha,3) + 1.7594e-6*pow(alpha,4));
71
72 delta_Cm_de = (-0.014928 - 0.0037783*alpha + 0.00039086*pow(de,2)
73 - 1.1304e-5*pow(de,3) - 1.8439e-6*pow(de,4));
74
75 delta_Cm = delta_Cm_a + delta_Cm_de;
76 KCm_alpha = delta_Cm_a/eta_ref_wing;
77 KCm_de = delta_Cm_de/eta_ref_tail;
78 delta_Cm = (0.75*eta_wing + 0.25*eta_tail)*KCm_alpha + (eta_tail)*KCm_de;
79
80 // hinge moment
81 if (alpha < 13)
82 {
83 delta_Ch_a = (-0.0012862 - 0.00022705*alpha + 1.5911e-5*pow(alpha,2)
84 + 5.4536e-7*pow(alpha,3));
85 }
86 else
87 {
88 delta_Ch_a = 0;
89 }
90 delta_Ch_e = -0.0011851 - 0.00049924*de;
91 delta_Ch = -(delta_Ch_a + delta_Ch_e);
92 KCh = -delta_Ch/eta_ref_tail;
93 delta_Ch = eta_tail*KCh;
94
95 // rolling moment
96 CL_diff = (eta_wing_left - eta_wing_right)*KCL;
97 delta_Cl = CL_diff/8.; // 10-23-02 Previously 4
98
99 //yawing moment
100 CD_diff = (eta_wing_left - eta_wing_right)*KCD;
101 delta_Cn = CD_diff/8.;
102
103 }
104
add_ice_effects()105 void add_ice_effects()
106 {
107 CD_clean = -1*CX*Cos_alpha*Cos_beta - CY*Sin_beta - CZ*Sin_alpha*Cos_beta;
108 CY_clean = -1*CX*Cos_alpha*Sin_beta + CY*Cos_beta - CZ*Sin_alpha*Sin_beta;
109 CL_clean = CX*Sin_alpha - CZ*Cos_alpha;
110 Cm_clean = Cm;
111 Cl_clean = Cl;
112 Cn_clean = Cn;
113 Ch_clean = Ch;
114
115 CD_iced = CD_clean + delta_CD;
116 CY_iced = CY_clean;
117 CL_iced = CL_clean + delta_CL;
118 Cm_iced = Cm_clean + delta_Cm;
119 Cl_iced = Cl_clean + delta_Cl;
120 Cn_iced = Cn_clean + delta_Cn;
121 //Ch_iced = Ch_clean + delta_Ch;
122
123 CD = CD_iced;
124 CY = CY_iced;
125 CL = CL_iced;
126 Cm = Cm_iced;
127 Cl = Cl_iced;
128 Cn = Cn_iced;
129 //Ch = Ch_iced;
130
131 CX = -1*CD*Cos_alpha*Cos_beta - CY*Cos_alpha*Sin_beta + CL*Sin_alpha;
132 CY = -1*CD*Sin_beta + CY*Cos_beta;
133 CZ = -1*CD*Sin_alpha*Cos_beta - CY*Sin_alpha*Sin_beta - CL*Cos_alpha;
134
135 }
136