付録B 主翼平面形最適化プログラム

/* 主翼平面形最適化プログラム(たんせい七號設計用) */

#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;
}

前の章 目次に戻る 次の章