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

#include
#include
#include
#include

#define JOINTMAX 10
#define LEN(r, t)   (Joint_pos[t] - Joint_pos[r])

#ifndef MAXDOUBLE
#include
#define MAXDOUBLE DBL_MAX
#endif

#ifndef PI
#define PI 3.141592
#endif

double Span;
double Total_area;
int Joint_count;
double Joint_pos[JOINTMAX];
double Chord_max;
double Chord_min;
int Taper_start;
double Step;
double Chord[JOINTMAX];
double Chord_opt[JOINTMAX];
double Air_density;
double Velocity;
double Lift_coefficient;
double Weight;
double Mean_chord_max;
double Mean_chord_min;
double Z_min;


void read_parameters()
{
 int i;

 printf("Wing Span [m] = ");
 scanf("%lf", &Span);

 printf("Wing Area [m^2] = ");
 scanf("%lf", &Total_area);

 printf("Number of Joints = ");
 scanf("%ld", &Joint_count);
 Joint_count += 2;

 Joint_pos[0] = 0;
 for(i = 1; i < Joint_count-1; i++) {
   printf("Joint Position(%u) [m] = ", i);
   scanf("%lf", &Joint_pos[i]);
 }
 Joint_pos[Joint_count-1] = Span / 2;

 printf("Chord maximum [m] = ");
 scanf("%lf", &Chord_max);

 printf("Chord minimum [m] = ");
 scanf("%lf", &Chord_min);

 printf("Taper start (1-%u) = ", Joint_count - 2);
 scanf("%ld", &Taper_start);

 printf("Calculation step [m] = ");
 scanf("%lf", &Step);

 printf("Mean chord max [m] = ");
 scanf("%lf", &Mean_chord_max);

 printf("Mean chord min [m] = ");
 scanf("%lf", &Mean_chord_min);

 printf("Weight [kg] = ");
 scanf("%lf", &Weight);

 printf("Velocity [m/s] = ");
 scanf("%lf", &Velocity);

 Air_density = 0.119;
 Lift_coefficient = 1.1;
 Z_min = MAXDOUBLE;

 printf("n");
}

double get_chord(double pos)
{
 int i;

 for(i = 0; Joint_pos[i+1] < pos; i++)
   ;
 pos -= Joint_pos[i];
 return Chord[i] + (Chord[i+1]-Chord[i]) * pos / LEN(i, i+1);
}

double integrate(double func(), double start, double end, int n)
{
 double h, s, x;

 h = (end - start) / (2 * n);
 x = start + h;
 s = func(start);
 while(--n) {
   s += 4 * func(x) + 2 * func(x + h);
   x += 2 * h;
 }
 s += func(end);

 return h * s / 3;
}

double circulation(double y)
{
 return (1 + 4 * Total_area * sqrt(1 - pow(2 * y / Span,2))
  / (PI * Span * get_chord(y))) * get_chord(y) * Velocity *     Lift_coefficient / 4;
}

double circulation_opt(double y)
{
 return 4 * Weight / (Air_density * Velocity * PI * Span) *  sqrt(1-pow(2 * y / Span, 2));
}

double evaluate_func(double y)
{
 return pow(circulation(y) - circulation_opt(y), 2);
}

double mean_chord_func(double y)
{
 return get_chord(y) * get_chord(y);
}

void evaluate()
{
 double z, mean_chord;
 int i;

/*  printf("r");
 for(i=0; i < Joint_count; i++)
   printf("%.3f ", Chord[i]);
*/
 z = integrate(evaluate_func, 0, Span / 2, 1000);
 mean_chord = 2 * integrate(mean_chord_func, 0, Span/2, 1000) /  Total_area;

 if (z < Z_min && Mean_chord_min < mean_chord && mean_chord <  Mean_chord_max) {
   memcpy(Chord_opt, Chord, sizeof(Chord));
   Z_min = z;
 }
}

void optimize(int pos, double area)
{
 double cr, ct, ct_max, ct_min;

 cr = Chord[pos];

 if((Span / 2 -Joint_pos[pos]) * cr < area)
   return;

 if(pos == 0) ct_max = cr;
 else ct_max = cr + (Chord[pos]-Chord[pos-1]) * LEN(pos,  pos+1)/LEN(pos-1, pos);

 if(pos == Joint_count - 2) {
   ct = 2 * area / LEN(pos, pos+1) - cr;
   if(ct < Chord_min || ct > ct_max) return;
   Chord[pos+1] = ct;
   evaluate();
   return;
 }

 if(pos < Taper_start) ct_min = cr;
 else ct_min = Chord_min;

 for(ct = ct_max; ct >= ct_min; ct -= Step) {
   Chord[pos+1] = ct;
   optimize(pos + 1, area - (cr + ct) * (Joint_pos[pos+1] -  Joint_pos[pos]) / 2);
 }
}

void show_result()
{
 int i;
 double pos;

 memcpy(Chord, Chord_opt, sizeof(Chord));

 printf("n------------------------------n");
 for(i = 0; i < Joint_count; i++)
   printf("Chord (%u) [mm] = %.3fn", i, Chord_opt[i]);
 printf("------------------------------n");
 printf("mean chord [mm] = %.3fn",
 2 * integrate(mean_chord_func, 0, Span/2, 1000) /  Total_area );
 printf("Z = %.3fn", Z_min);
 printf("------------------------------n");
 printf("Circulation distributionn");
 for(pos = 0; pos <= Span/2; pos += 0.01)
   printf("%.2f %.3fn", pos, circulation(pos));
 printf("------------------------------n");
 printf("Cl distributionn");
 for(pos = 0; pos <= Span/2; pos += 0.01)
   printf("%.2f %.3fn", pos, 2 * circulation(pos) /  (get_chord(pos) * Velocity));
}


int main()
{
 double c;

 read_parameters();
 for(c = Total_area / Span; c <= Chord_max; c += Step) {
   Chord[0] = c;
   optimize(0, Total_area / 2);
 }
 show_result();

 return 0;
}

前の章 目次に戻る 次の章