-
Notifications
You must be signed in to change notification settings - Fork 9
/
longtrim.cpp
256 lines (224 loc) · 9.83 KB
/
longtrim.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
void LongitudinalTrim(GENERAL,PANEL *,DVE *,int,double *&,double &,\
double &,FILE *,double ***);
void LongitudinalTrim(GENERAL info,PANEL *panelPtr,DVE *surfaceDVEPtr,int HTpanel,\
double *&cn,double &CL,double &CDi,\
FILE *MomSol,double ***camberPtr)
{
// This program finds longitudinal trim solutions for twin wing configurations
//
// Wing 1 is the main wing, Wing 2 is the horizontal tail.
//
//ATTENTION: NO SPANWISE VELOCITY VARIATION !!!
//If you don't know what this means, don't worry and move on.
//
// INPUT
// info general information
// panelPtr panel information
// HTpanel index of first panel of HT
//
//output
// surfaceDVEPtr geometry of surface elements
// cn normal force coefficient of each span element
// cd drag force coefficient of each span element
// CL total lift
// CDi total induced drag
//
//the routine also saves the trim information to TrimSol.txt
//furthermore saved is the spanwise lift distribution and spanwise
//drag distribution in SpanAOA<info.alpha>.txt where <info.alpha>
//holds the numerical value of the current angle of attack
int i,j,k,l,m,n=0; //loop counter
double epsilonHT; //HT incident correction [rad]
double CM_resid; //residual moment
double eps1; //the HT incidence angle of one iteration previous
double CM_old; //saving pitching moment of previous iteration step
double CLht,CLhti; //total and induced lift of HT
double *cl,*cy; //section lift and side force coefficients
double *S; //area of a spanwise strip (sum of areas of DVEs of one span location)
double tempS;
int tempI;
double **N_force; //surface DVE's normal forces/density
//[0]: free stream lift, [1]: induced lift,
//[2]: free stream side, [3]: induced side force/density
//[4]: free str. normal, [5]: ind. normal frc/density
//[6,7,8]: eN_x, eN_y, eN_z in global ref. frame
double *D_force; //drag forces/density along span
double **Span_force; //x,y,z aerodynamic force/density in wind-axis system
double *cd; //section ind. drag coefficient
FILE *spaninfo; //output file for spanwise information
char filename[133]; //file path and name for spanwise information
//allocating memory
ALLOC1D(&cl,info.nospanelement); //section lift coefficient
ALLOC1D(&cy,info.nospanelement); //section side force coefficient
ALLOC1D(&S,info.nospanelement); //section area
ALLOC1D(&cd,info.nospanelement); //section ind. drag coefficient
ALLOC2D(&N_force,info.noelement,9); //surface DVE normal forces
ALLOC1D(&D_force,info.nospanelement);//Drag force per span element
ALLOC2D(&Span_force,info.nospanelement,3);//Span force vector per span element
//initial HT incident correction
epsilonHT = 0;//panelPtr[i].eps1; //[rad]
eps1 = 0; //the HT incidence angle of one iteration previous
//===========================================//
//pitching moment computation - Step 0
//===========================================//
//computed the residual moment coefficient of wing and tail
CM_resid = PitchingMoment\
(info,panelPtr,surfacePtr,info.cmac,epsilonHT,\
HTpanel,info.RefPt,CLht,CLhti,\
N_force,D_force,Span_force,CL,CDi,camberPtr);
//Subroutine in PitchMoment.cpp
//adding zero lift of wing only if camber is turned off
if(~info.flagCAMBER){CM_resid += info.CMoWing;}
if(info.trim==1) //longitudinal trim routine
{
if(CM_resid > 0) epsilonHT = 0.035; //deflect HT trailing edge 2deg down
else epsilonHT =-0.035; //deflect HT trailing edge 2deg up
//========================================//
// pitching moment iteration loop
//========================================//
i=0; //initializing loop counter
do
{
//========================================//
// pitching moment computation - Step i
//========================================//
//saving pitching moment of previous iteration step
CM_old = CM_resid;
//computed the residual moment coefficient of wing and tail
CM_resid = PitchingMoment\
(info,panelPtr,surfacePtr,info.cmac,epsilonHT,\
HTpanel,info.RefPt,CLht,CLhti,\
N_force,D_force,Span_force,CL,CDi,camberPtr);
//Subroutine in PitchMoment.cpp
//adding zero lift of wing
CM_resid += info.CMoWing;
//computing new HT incidence angle
tempS = epsilonHT - CM_resid*(epsilonHT-eps1)/(CM_resid-CM_old);
eps1 = epsilonHT; epsilonHT = tempS; //reassigning HT angles
i++; //incrementing loop counter
}while((i<20) && (CM_resid*CM_resid>.000025));
//========================================//
// pitching moment iteration loop END
//========================================*///
} //end of longitudinal trim routine
//===================================================================//
// Compute cn and cd (spanwise values)
//===================================================================//
i=0; //intitaliing span index counter
m=0; //index of leading edge DVE
for(k=0;k<info.nopanel;k++) //loop over panels
{
for(l=0;l<panelPtr[k].n;l++) //loop over span of panel k
{
//working each spanwise strip
//initializing
cn[i]=0;
cl[i]=0;
cy[i]=0;
S[i] = 0; //area of current spanwise strip
//adding up chordal values of one spanwise location (indexed i)
for(m=0;m<panelPtr[k].m;m++)
{
j=n+m*panelPtr[k].n; //counting index along chord
//
cn[i] += (N_force[j][4]+N_force[j][5]); //adding forces/density along chord
cl[i] += (N_force[j][0]+N_force[j][1]);
cy[i] += (N_force[j][2]+N_force[j][3]);
S[i] += surfaceDVEPtr[j].S; //adding DVE areas along chord
}
//Nondimensionalizing values using summed areas and velocity at LE of chordal row of DVEs
tempI = j-panelPtr[k].n*(panelPtr[k].m-1); // index of DVE at leading edge
tempS = 2/(S[i]*dot(surfaceDVEPtr[tempI].u,surfaceDVEPtr[tempI].u));
cd[i] = D_force[i]*tempS;
cn[i] *= tempS;
cl[i] *= tempS;
cy[i] *= tempS;
i++; //next span index
n++; //index of next leading edge DVE
}
n += panelPtr[k].n*(panelPtr[k].m-1); //index of next LE DVE of next panel
}
for(i=0;i<info.nospanelement;i++)
{
printf("%d cl %lf cy %lf cn %lf cd %lf in longtrim\n",i,cl[i],cy[i],cn[i],cd[i]);
}
for(i=0;i<info.nospanelement;i++)
{
printf("%d cfz %lf cfy %lf cfx %lf in longtrim\n",i,Cf[i][2],Cf[i][1],Cf[i][0]);
}
//===================================================================//
// save spanwise information, lift and drag distribution
//===================================================================//
//creates file name AOA#.##.txt ## is angle of attack
sprintf(filename,"%s%s%.2lf%s",OUTPUT_PATH,"AOA",info.alpha*RtD,".txt");
//creates file in subdirectory output
spaninfo = fopen(filename, "w");
//write header
fprintf(spaninfo,"This file contains spanwise information at");
fprintf(spaninfo," angle of %.2lf\n",info.alpha*RtD);
fprintf(spaninfo," The total lift (uncorrected for stall) and drag coefficient are ");
fprintf(spaninfo," CL = %6.6lf (uncorrected) CDi = %12.8lf ",CL,CDi);
fprintf(spaninfo," The horizontal tail is panel %d through %d\n"\
,info.wing1[1],info.wing2[1]);
//writes header for information on surface elements
fprintf(spaninfo,"%6s %16s %16s %16s %16s %16s %16s %16s",\
"index","xo","yo","zo","cn","cl","cy","cd");
fprintf(spaninfo," %16s %16s %16s %16s %16s %16s %16s %16s %16s %16s\t#\n",\
"A_TE","B_TE","C_TE","S","span","chord","nu","epsilon","psi","phiLE");
// for(i=0;i<info.nospanelement;i++)
i=0; //intitaliing span index counter
m=0; //index of leading edge DVE
for(k=0;k<info.nopanel;k++) //loop over panels
{
for(l=0;l<panelPtr[k].n;l++) //loop over span of panel k
{
if(i==info.wing1[1])
fprintf(spaninfo,"HT\n"); //separates HT data
j=m+(panelPtr[k].n*(panelPtr[k].m-1));
//surface element index
fprintf(spaninfo,"%6d",i);
//coord. of ref point
fprintf(spaninfo," %16.12lf %16.12lf %16.12lf",\
surfacePtr[m].xo[0],surfacePtr[m].xo[1],\
surfacePtr[m].xo[2]);
//normal, lift, side, drag force coefficients
fprintf(spaninfo," %16.12lf %16.12lf %16.12lf %16.12lf",\
cn[i],cl[i],cy[i],cd[i]);
//more info on element
fprintf(spaninfo," %16.12lf %16.12lf %16.12lf %16.12lf %16.12lf %16.12lf",\
surfacePtr[j].A,surfacePtr[j].B,surfacePtr[j].C,S[i],\
surfacePtr[m].eta*2,surfacePtr[m].xsi*2*panelPtr[k].m);
fprintf(spaninfo," %16.12lf %16.12lf %16.12lf",\
surfacePtr[m].nu*RtD,surfacePtr[m].epsilon*RtD,\
surfacePtr[m].psi*RtD);
fprintf(spaninfo," %16.12lf",\
(surfacePtr[m].phiLE*RtD));
fprintf(spaninfo,"\n");
i++; //next span index
m++; //index of next leading edge DVE
}
m += panelPtr[k].n*(panelPtr[k].m-1); //index of next LE DVE of next panel
}
fclose(spaninfo);
//===================================================================//
// DONE saving spanwise information, lift and drag distribution
//===================================================================//
//===================================================================//
//write trim results to TrimSol.txt
//===================================================================//
fprintf(MomSol,"%6.2lf %6.2lf ",info.alpha*RtD,epsilonHT*RtD);
fprintf(MomSol,"%6.3lf %12.8lf %7.4lf ",CL,CDi,CM_resid);
fprintf(MomSol,"%10.6lf %8.4lf\n",CLht,info.CMoWing);
fflush(MomSol);
FREE2D(&N_force,info.noelement,6);
FREE1D(&D_force,info.nospanelement);
FREE2D(&Span_force,info.nospanelement,3);
FREE1D(&cl,info.nospanelement); //section lift coefficient
FREE1D(&cy,info.nospanelement); //section side force coefficient
FREE1D(&S,info.nospanelement); //section area
FREE1D(&cd,info.nospanelement); //section side force coefficient
printf("\n");
}
//===================================================================//
//END of program
//===================================================================//