/* 主翼平面形最適化プログラム(たんせい七號設計用) */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#ifndef PI
#define PI 3.1415926535
#endif
#define TRUE 1
#define FALSE 0
#define EPS 1.0e-6
#define M 100
#define TH(n) (((n)+1) * PI / (M + 1))
#define MAX_Y_CNT 10
#define AIR_DENSITY 0.119
double Weight;
double Velocity;
double Span;
double Y[MAX_Y_CNT];
double Chord[MAX_Y_CNT];
double Wing_area;
double Aspect_ratio;
double Angle_of_attack;
double A0;
double Circulation_opt[M];
double Chord_opt[MAX_Y_CNT];
double Bb[M][M];
double Step;
double Chord_min;
int Y_cnt;
int Taper_start;
int Found = FALSE;
void read_parameters()
{
int i;
printf("機体重量[kg] = ");
scanf("%lf", &Weight);
printf("巡行速度[m/s] = ");
scanf("%lf", &Velocity);
printf("主翼面積[m^2] = ");
scanf("%lf", &Wing_area);
printf("コード長を変化させる位置の数 = ");
scanf("%d", &Y_cnt);
for(i=0; i < Y_cnt; i++) {
printf("コード長を変化させる位置(%d番目)[m] = ", i);
scanf("%lf", &Y[i]);
}
Span = 2 * Y[Y_cnt - 1];
Aspect_ratio = pow(Span, 2) / Wing_area;
printf("テーパーを開始する位置(0-%d) = ", Y_cnt - 1);
scanf("%d", &Taper_start);
printf("迎角(揚力=0となる角度から測る)[degree] = ");
scanf("%lf", &Angle_of_attack);
Angle_of_attack = Angle_of_attack / 180 * PI;
A0 = 2 * PI;
printf("最適化時にコード長を変化させるきざみ[m] = ");
scanf("%lf", &Step);
printf("コード長の最小値[m] = ");
scanf("%lf", &Chord_min);
}
void init_table()
{
int i, j;
for(i = 0; i < M; i++) {
for(j = 0; j < M; j++) {
if(i == j) Bb[i][j] = (M + 1) / (4 * sin(TH(i)));
else Bb[i][j] = (1 - pow(-1, i - j)) / (2 * (M + 1)) * sin(TH(j)) / pow(cos(TH(i)) - cos(TH(j)),2);
}
}
}
double get_chord(double pos)
{
int t;
for(t = 0; Y[t] < pos; t++);
return Chord[t - 1] + (pos - Y[t - 1]) * (Chord[t] - Chord[t - 1]) / (Y[t] - Y[t - 1]);
}
void calc_circulation(double *circulation)
{
int i, j, again;
double b[M], t;
for(i = 0; i < M; i++)
b[i] = Bb[i][i] + 2 * Span / (A0 * get_chord(fabs(Span / 2 * cos(TH(i)))));
for(i = 0; i < M; i++)
circulation[i] = 2 / (PI * Aspect_ratio) * sqrt(1 - pow(cos(TH(i)),2));
do {
again = FALSE;
for(i = 0; i < M; i++) {
t = Angle_of_attack - Bb[i][i] * circulation[i];
for(j = 0; j < M; j++)
t += Bb[i][j] * circulation[j];
t /= b[i];
if(fabs(t - circulation[i]) > EPS) again = TRUE;
circulation[i] = t;
}
} while(again);
}
double circulation_goal(double y)
{
double ret;
return ret = 4 * Weight / (AIR_DENSITY * Velocity * PI * Span) * sqrt(1 - pow(2 * y / Span, 2));
}
void evaluate()
{
static double z_min = 0;
static int firsttime = TRUE;
double circulation[M], z, grad, t;
int i;
if(Chord[Y_cnt - 1] < Chord_min) return;
for(i = 0, grad = 0; i < Y_cnt - 1; i++) {
t = (Chord[i] - Chord[i + 1]) / (Y[i + 1] - Y[i]);
if(t < grad) return;
grad = t;
}
calc_circulation(circulation);
for(z = 0, i =0; i < M; i++)
z += pow(circulation[i] - circulation_goal(Span / 2 * cos(TH(i))), 2) * sin(TH(i));
if(firsttime || z < z_min) {
firsttime = FALSE;
z_min = z;
memcpy(Chord_opt, Chord, sizeof(Chord));
memcpy(Circulation_opt, circulation, sizeof(circulation));
Found = TRUE;
}
}
void optimize(int pos, double area)
{
double lt, lp, cp;
lt = Span / 2 - Y[pos];
lp = Y[pos + 1] - Y[pos];
if(Chord[pos] * lt < area || Chord[pos] * lt / 2 > area) return;
if(pos + 1 == Y_cnt - 1) {
Chord[pos + 1] = 2 * area / lp - Chord[pos];
evaluate();
} else {
if(pos < Taper_start) {
Chord[pos + 1] = Chord[pos];
optimize(pos + 1, area - Chord[pos] * lp);
} else {
for(cp = Chord[pos]; cp >= 0; cp -= Step) {
Chord[pos + 1] = cp;
optimize(pos + 1, area - (Chord[pos] + cp) * lp / 2);
}
}
}
}
void print_result()
{
int i;
if(Found) {
printf("nコード長n");
printf("Y[m]tChord[m]n");
for(i = 0; i < Y_cnt; i++)
printf("%ft%fn", Y[i], Chord_opt[i]);
printf("n循環分布(スパン×機速で無次元化した値)n");
printf("位置[m]t循環n");
for(i = 0; i < M; i++)
printf("%ft%fn", cos(TH(i)) * Span / 2, Circulation_opt[i]);
} else {
printf("最適化できませんでしたn");
}
}
void optimize_root()
{
double c0;
for(c0 = 2 * Wing_area / Span; c0 > 0; c0 -= Step) {
Chord[0] = c0;
optimize(0, Wing_area / 2);
}
}
int main()
{
read_parameters();
init_table();
optimize_root();
print_result();
return 0;
}
前の章 目次に戻る 次の章