/*********************************************************************/
/* Calcul des profils aerodynamiques par la méthode des panneaux, */
/* avec tourbillons et par formulation des lignes de courant */
/* --------------------- session automne 2000 ---------------------- */
/* */
/* Ecole polytechnique de Montreal */
/* */
/* Conception et realisation: */
/* */
/* BOLOVIS Angela */
/* FOISSAC Émilie */
/* LEROUYER Benjamin */
/* MEUNIER Arnaud */
/* */
/*********************************************************************/
/* inclusion dans le programme des fonctions utilisees */
#include <stdio.h>
#include <string.h>
#include <dos.h>
#include <math.h>
#include <stdlib.h>
#include <conio.h>
/* Définition de la constante qui nous permettra de convertir des degres en radians*/
const double Pi = 3.141592653589793238462643383279502884197169;
const double deg_rad = (3.141592653589793238462643383279502884197169/180);
const double pi_fois_2 = 3.141592653589793238462643383279502884197169*2;
/* definition du nombre de panneaux entres par l utilisateur comme variable globale */
int N;
/*********************************************************************/
/* Les references aux formules du livre de Ion Paraschiviu de 1998, */
/* aerodynamique subsonique, seront notees [4.130] par exemple */
/*********************************************************************/
/**************************************************************************/
/* Variables utilisees souvent dans le programme: */
/* N = nombre de panneaux demandes par l'utilisateur */
/* A noter que le nombre de panneaux = nombre de noeud car le */
/* profil est ferme */
/* A = matrice carree des coefficients d influence, de taille N+1 */
/* b = vecteur de taille N+1 (vecteur des vitesses avant la résolution */
/* du systeme lineaire (Gauss), puis vecteur des tourbillons ensuite */
/* Xi et Yi sont les tableaux contenant la position des noeuds i */
/* X et Y sont les tableaux sur la position des points de controles */
/* a noter que les X, Y, Xi et Yi sont alloues avec une taille de */
/* N+1, afin que le panneau n+1 ai les memes coordonnees */
/* que le panneau 1, pour simplifier les algorithmes. */
/* r est le rayon du bord d'attaque en relatif */
/* d est la position relative de la cambrure max */
/* f est la cambrure relative max */
/* t est l epaisseur relative max */
/* */
/* Lors de l'utilisation d'une matrice, nous nous servirons des indices */
/* i et j tels que "tableau[i][j]", avec i la ligne et j la colonne. */
/**************************************************************************/
/*********************************************************************/
/* Saisie des donnees necessaires au calcul et interface utilisateur */
/*********************************************************************/
void menu (int *naca, double *v_inf, double *alpha)
{
char reponse;
reponse = 'n';
while (reponse != 'o' && reponse != 'O')
{
/* verification de la validite du profil NACA rentre par utilisateur
on ne prend que les profils a 4 chiffres */
*naca = 10000;
while (*naca > 9999)
{
printf ("\nQuel est le numero NACA du profil 4 chiffres a etudier ?\n");
scanf ("%d",&*naca);
if (*naca > 9999) printf ("\n Ceci ne correspond pas à un profil a 4 chiffres");
}
do{
printf ("\nCombien souhaitez vous de panneaux ? (60 panneaux permettent une assez bonne precision) : ");
scanf ("%i",&N);
}while (N<3);
printf ("\nQuelle est la vitesse a l'infini (em m/s)? : ");
scanf ("%lf",&*v_inf);
printf ("\nQuel est l'angle d'incidence du profil (en degres) ? : ");
scanf ("%lf",&*alpha);
if (*naca < 100) printf ("\n\nVous avez demande l'etude d'un profil NACA symetrique 00%d\n", *naca);
else printf ("\n\nVous avez demande l'etude d'un profil NACA cambre %d\n", *naca);
printf ("\nCes donnees vous conviennent elles ? (o/n)\n");
fflush(stdin);
scanf ("%c",&reponse);
}
fflush(stdin);
}/* fin de la fonction menu */
/*********************************************************************/
/* Saisie des donnees necessaires aux calculs apres premier passage */
/*********************************************************************/
int menu2 (int naca, double *alpha, double *alpha2)
{
int reponse;
clrscr();
printf ("\n\nLe calcul avec les valeurs precedemment entrees s'est bien deroule. Voulez-Vous? ");
printf ("\n\n 1- Quitter le programme ?");
printf ("\n 2- Relancer les traitements avec un nombre de panneaux different ?");
printf ("\n 3- Modifier le nombre de panneaux et l'angle d'incidence ?");
printf ("\n 4- Modifier le nombre de panneaux et calculer la courbe Cl\n");
printf (" en fonction de alpha ?");
printf ("\n\n Entrez votre choix :");
fflush(stdin);
scanf ("%d",&reponse);
fflush(stdin);
if (reponse == 1) return (reponse);
printf ("\nCombien souhaitez vous de panneaux ? (60 panneaux permettent une assez bonne precision) : ");
scanf ("%i",&N);
if (reponse == 2) return (reponse);
if (reponse == 3)
{
printf ("\nQuel est le nouvel angle d'incidence du profil (en degres) ? : ");
scanf ("%lf",&*alpha);
return (reponse);
}
printf ("\nQuel est l'angle d'incidence du profil au depart de la courbe Cl(alpha) (en degres) ? : ");
scanf ("%lf",&*alpha);
do
{
printf ("\nQuel est l'angle d'incidence du profil a la fin de la courbe Cl(alpha) (en degres) ? : ");
scanf ("%lf",&*alpha2);
}while (*alpha2<*alpha);
return (reponse);
}/* fin de la fonction menu */
/**************************************************************************/
/* Calcul des caractéristiques du profil a partir du NACA rentre */
/* Le premier chiffre du NACA donne la cambrure relative max (f) */
/* Le deuxième chiffre donne la position relative de la cambrure max (d) */
/* Les deux derniers indiquent l épaisseur relative max (t) */
/**************************************************************************/
void caracteristiques (int naca, char *sym, double *d, double *f, double *r, double *t)
{
int d_bis, f_bis, t_bis;
/* les variables bis nous permettent de travailler avec le NACA de type entier */
f_bis = (naca / 1000);
d_bis = ((naca % 1000) /100);
t_bis = (naca % 100);
/* On remplis les valeurs caracteristiques relatives du profil */
*d = d_bis;
*d = *d / 10;
*f = f_bis;
*f = *f /100;
*t = t_bis;
*t = *t / 100;
*r = 1.1019*(*t)*(*t);
/* test de la symétrie du profil*/
if (naca / 100 == 0) *sym = 's';
else *sym = 'c';
}
/*******************************************************************/
/* Calcul des coordonnées des points de discrétisation */
/* Cette fonction est dépendante de l équation du profil étudie, */
/* ici le NACA. On ne parcourra qu'une fois la corde, en calculant*/
/* pour chaque position de x les ordonnées sur l extrados et sur */
/* l intrados. */
/* les calculs seront différents si N est pair ou impair. */
/* Le dernier nœud de l'intrados (cote bord d'attaque) sera */
/* toujours N/2, en prenant la valeur inférieur de ce résultat si */
/* N est impair. Si N est pair, il y aura un nœud situe juste sur */
/* le bord d attaque. Les nœuds de l'extrados ont la même abscisse*/
/* que les nœuds correspondant de l'intrados */
/*******************************************************************/
void discretisationNACA (char sym, double f, double d, double t, double * Xi, double *Yi)
{
int i, nb=N+1;
double yc, teta;
/* le premier noeud est toujours a la position (1,0) */
Xi[0] = 1;
Xi[N] = 1;
Yi[0] = 0;
Yi[N] = 0;
/* si N est pair le dernier noeud de l'intrados est N/2
et le noeud (N/2)+1 est sur le bord d'attaque, donc a la position (0,0) */
if (N%2==0)
{
Xi[(N/2)] = 0;
Yi[(N/2)] = 0;
nb--;
}
/* sinon si N est impair, le dernier point de l'intrados est ((N-1)/2)+1 et il n'y a pas de
noeud a la position (0,0) */
/* ---- Si le profil est symetrique on ne calcule que la moitie des noeuds
les autres etant trouves par symetrie d axe Ox */
if (sym == 's')
{
/* on se deplace du noeud 2 au dernier noeud de l'intrados (varie suivant que N est paire ou impaire)*/
for (i=1; i<(nb/2); i++)
{
Xi[i] = 0.5*(1 + cos (2*Pi*i/N));
/* equation du profil etudie (dans notre cas un profil NACA) */
Yi[i] = -(t*5)*(0.2969*sqrt(Xi[i]) - 0.126*Xi[i] - 0.3537*Xi[i]*Xi[i] + 0.2843*Xi[i]*Xi[i]*Xi[i] - 0.1015*Xi[i]*Xi[i]*Xi[i]*Xi[i]);
/* Le noeud N est a la meme abscisse que le noeud 2, N-1 que 3 ... donc
le noeud (N+2)-i et le noeud i ont meme abscisse */
Xi[(N-i)] = Xi[i];
Yi[(N-i)] = -Yi[i];
}
}
/* ---- si le profil n est pas symetrique */
else
{
/* on se deplace du noeud 2 au dernier noeud de l'intrados (varie suivant que N est paire ou impaire)*/
for (i=1; i<(nb/2); i++)
{
/* calcul du profil symetrique equivalent */
Xi[i] = 0.5*(1 + cos (2*Pi*i/N));
/* Yi[i] = (t/0.2)*(0.2969*sqrt(Xi[i]) + Xi[i]*(-0.126 + Xi[i]*(-0.3537 + Xi[i]*(0.2843 - 0.1015*Xi[i])))); */
Yi[i] = (t*5)*(0.2969*sqrt(Xi[i]) - 0.126*Xi[i] - 0.3537*Xi[i]*Xi[i] + 0.2843*Xi[i]*Xi[i]*Xi[i] - 0.1015*Xi[i]*Xi[i]*Xi[i]*Xi[i]);
/* calcul de l ordonnee de la ligne de cambrure moyenne */
if ((0 <= Xi[i]) && (Xi[i] <= d))
{
yc = (f*Xi[i]*(2*d-Xi[i]))/(d*d);
teta = deg_rad*(2*f*(d-Xi[i]))/(d*d);
}
else
{
yc = (f*(1-Xi[i])*(1+Xi[i]-2*d))/((1-d)*(1-d));
teta = deg_rad*(2*f*(d-Xi[i]))/((1-d)*(1-d));
}
/* coordonnees sur le profil cambre a partir de la ligne de cambrure moyenne*/
/* points sur l'extrados */
Xi[(N-i)] = Xi[i] - Yi[i] * sin(teta);
Yi[(N-i)] = yc + Yi[i] * cos(teta);
/* points sur l'intrados */
Xi[i] = Xi[i] + Yi[i] * sin(teta);
Yi[i] = yc - Yi[i] * cos(teta);
}
}
}/* fin de la fonction discretisation */
/*******************************************************************/
/* Allocation dynamique d'une matrice */
/* la matrice passee en parametre est cree en memoire, avec une */
/* taille de nblig lignes et nbcol colonnes. L adresse du premier */
/* pointeur de la matrice est retourne */
/*******************************************************************/
double ** alloc_mat (double **mat, int nblig, int nbcol)
{
int i;
/* allocation d`un tableau de nblig pointeurs dont chaque element va pointer
sur un tableau de nbcol elements, ce tableau representant une ligne */
if ((mat=(double**)malloc((nblig)*sizeof(double *)))==NULL)
{
printf ("allocation ligne impossible");
getch();
exit(0);
}
else
{
for (i=0;i<=(nblig-1);i++)
{
if ((mat[i]=(double*)malloc((nbcol)*sizeof(double)))==NULL)
{
printf("allocation colonne impossible");
getch();
exit(0);
}
}
}
return mat;
}/* fin de alloc_mat */
/*******************************************************************************/
/* Allocation dynamique d'un vecteur */
/*******************************************************************************/
double* alloc_vec(double * vect, int nblig)
{
if ((vect = (double*)malloc((nblig)*sizeof(double)))==NULL)
{
printf("allocation impossible");
getch();
exit(1);
}
return vect;
}/* fin de alloc_vec */
/******************************************************************************/
/* Remplissage de la matrice A des coefficients */
/******************************************************************************/
void initmatrice (double **A, double *Xi, double *Yi, double *X, double *Y, double *l)
{
int i,j;
double xetoil, yetoil, betapj, **r, *costeta, *sinteta, dxi, dyi, dxijp, dyijp;
/* Allocation des vecteurs et matrice necessaires aux calculs */
sinteta = alloc_vec (sinteta,N); /* inclinaison des panneaux (sin teta j) */
costeta = alloc_vec (costeta,N); /* inclinaison des panneaux (cos teta j) */
r = alloc_mat (r,N,N+1); /* distance du noeud j au milieu du panneau i */
/* calcul des x et y moyens (coordonnees des points de controle) [4.130] */
for (i=0;i<N;i++)
{
X[i] = (Xi[i] + Xi[i+1])/2;
Y[i] = (Yi[i] + Yi[i+1])/2;
}
/* le point de controle du panneau N+1 est celui du panneau 1 */
X[N] = X[0];
Y[N] = Y[0];
for (i=0;i<N;i++)
{
/* calcul des rij [4.157]*/
for (j=0;j<N+1;j++)
{
r[i][j] = sqrt( ((X[i]-Xi[j])*(X[i]-Xi[j])) + ((Y[i]-Yi[j])*(Y[i]-Yi[j])) );
}
/* composantes du panneau i */
dxi = Xi[i+1]-Xi[i];
dyi = Yi[i+1]-Yi[i];
/* calcul de l, le tableau vecteur des longueurs des panneaux */
l[i] = sqrt( (dxi*dxi) + (dyi*dyi) );
/* calcul des teta j */
costeta[i] = (dxi/l[i]);
sinteta[i] = (dyi/l[i]);
}
/* Calcul des coefficients Aij de la matrice A */
for (i=0;i<N;i++)
{
for (j=0;j<N;j++)
{
/* Remplissage de la diagonale [4.178]*/
if (i==j) A[i][j] = l[j] * (log(l[j]/2)-1) / pi_fois_2;
else
/* application de l equation [4.175] */
{
dxi = X[i]-Xi[j];
dyi = Y[i]-Yi[j];
dxijp = X[i]-Xi[j+1];
dyijp = Y[i]-Yi[j+1];
/* Definition de xetoil et yetoil [4.176 et 4.177]*/
xetoil = (dxi*costeta[j]) + (dyi*sinteta[j]);
yetoil = (dyi*costeta[j]) - (dxi*sinteta[j]);
/* Definition de beta p,j , avec p=i , j=j donc Xp=Xi[p],Yp=Yi[i] :coordonnées du centre du panneau P. */
/* betapj = angle du panneau j vu du milieu du panneau i [4.158]*/
betapj = atan2( ((dxi*dyijp) - (dyi*dxijp)),((dyi*dyijp) + (dxi*dxijp)));
/* Initialisation des termes non diagonaux [4.175] */
// A[i][j] = ( (l[j]*(log(r[i][j+1])-1)) + (xetoil*log(r[i][j]/r[i][j+1])) + (yetoil*betapj) )/pi_fois_2;
A[i][j] = l[j] * ( log (r[i][j+1]) - 1 );
A[i][j] = A[i][j] + ( xetoil * log ( r[i][j] / r[i][j+1] ) );
A[i][j] = A[i][j] + (yetoil*betapj);
A[i][j] = A[i][j] / pi_fois_2;
}
/* initialisation des dernieres colonne et ligne de la matrice A */
A[N][j] = 0;
A[j][N] = -1;
}
}
/* initialisation des dernieres colonne et ligne de la matrice A */
A[N][N-1] = 1;
A[N][0] = 1;
A[N][N] = -1;
/* desallocations des tableaux */
free(costeta);
free(sinteta);
free(r);
} /* fin de initmatrice */
/*************************************************************************/
/* Remplissage du vecteur b des vitesses */
/*************************************************************************/
void init_b (double *b, double *v_inf, double alpha, double *X, double *Y)
{
int i;
/* N premières valeurs du vecteur b [4.182] */
for (i=0;i<N;i++)
{
b[i] = -*v_inf * (Y[i]*cos(deg_rad*alpha) - X[i]*sin(deg_rad*alpha));
}
/* N+1 eme valeur du vecteur b [4.183] */
b[N] = 0;
} /* fin de init_b */
/****************************************************************************/
/* Resolution d'un systeme lineaire par la methode du pivot de Gauss */
/* Cette fonction resout le systeme suivant: [A].{x}={b} avec la matrice A et
le vecteur b connus. Cette fonction nous renvoie le vecteur x dans la matrice
b passee en parametre */
/****************************************************************************/
int gauss_pivmax (double **a,double *b)
{
int i,j,k,l,err;
double max,pivot,coef,s;
double **t;
t = alloc_mat (t,N+1,N+2);
for (i=0;i<N+1;i++)
{
for (j=0;j<N+1;j++)
{
t[i][j] = a[i][j];
t[i][N+1] = b[i];
}
}
err = 1;
k = 0;
while (err==1 && k<N)
{
max = fabs(t[k][k]);
l = k;
for (i=k+1;i<N+1;i++)
{
if (max < fabs(t[i][k]))
{
max = fabs(t[i][k]);
l = i;
}
}
if (max!=0)
{
if (l!=k)
{
for (j=k;j<N+2;j++)
{
pivot = t[k][j];
t[k][j] = t[l][j];
t[l][j] = pivot;
}
}
pivot = t[k][k];
for (i=k+1;i<N+1;i++)
{
coef = t[i][k]/pivot;
for (j=k+1;j<N+2;j++) t[i][j] -= coef*t[k][j];
}
}
else err = 0;
k++;
} /* fin de boucle while */
/* Si le pivot de gauss est nul, on ne peut pas retourner de */
/* vecteur resultat. retour du code erreur 0 */
if (t[N][N]==0) err=0;
if (err==1)
{
b[N] = t[N][N+1]/t[N][N];
for (i=N-1;i>=0;i--)
{
s = t[i][N+1];
for (j=i+1;j<N+1;j++) s -= t[i][j]*b[j];
b[i] = s/t[i][i];
}
}
return(err);
} /* fin de la fonction du pivot de Gauss */
/*****************************************************************************/
/* Calcul du coefficients de pression */
/*****************************************************************************/
void C_pression (double *x, double *v_inf, double *Cp)
{
int i;
/* calcul des coefficients de pression [4.142] */
for (i=0;i<N;i++)
{
Cp[i] = 1 - (x[i] / *v_inf)*(x[i] / *v_inf);
}
}
/*****************************************************************************/
/* Calcul du coefficients de portance */
/*****************************************************************************/
void C_portance (double *x, double *v_inf, double *Cl, double *l)
{
int i;
Cl[N]=0;
for (i=0;i<N;i++)
{
/* calcul du coefficient de portance du panneau i , au point de controle [4.167]*/
Cl[i] = (2/(*v_inf))*(x[i]*l[i]) ;
/* calcul du coefficient de portance de l aile entiere*/
Cl[N] = Cl[N] + Cl[i];
}
}/* fin de la fonction de calcul du coefficient de portance */
/***********************************************************************/
/* Ecriture des résultats dans un fichier interface avec MAPLE, resultataero.txt
/* le fichier contiendra 18 enregistrements d'entete, suivi de N */
/* enregistrements contenant les informations du panneau correspondant */
/* i.e. au panneau i est associé le noeud i, le point de controle du */
/* panneau i ainsi que les coefficients de portance et de pression. */
/* avec i variant de 1 a N */
/* Un enregistrement contiendra les champs suivants:
double : abscisse point de controle noeud i (X)
double : ordonnee du point de controle noeud i (Y)
double : coefficient pression du noeud i (pression)
double : coefficient de portance du panneau i (portance)
double : abscisse du noeud i (Xi)
double : ordonnee du noeud i (Yi)
Un autre fichier sera genere, resultat.txt identique au premier, mais
sans l entete ni les separateurs de champs. Ce fichier sera rapatrie
dans Microsoft excel, ou mathlab, ou Maple, ou n importe quel
autre programme capable de generer des courbes.
Attention, sous excel, il faudra remplacer tous les separateurs
decimaux "." par les separateurs decimaux "," */
/**********************************************************************/
void ecrirefichier(double *X,double *Y,double *Xi,double *Yi,double *Cl,double *Cp,int naca, double *d, double *f, double *r, double alpha)
{
int i;
FILE *F=fopen("resultataero.txt","w");
FILE *F1=fopen("resultat.txt","w");
FILE *F2=fopen("resultat_X_Cp.txt","w");
FILE *F3=fopen("resultat_X_Y_Cp.txt","w");
/* si ouverture du fichier impossible*/
if ((F==NULL)||(F1==NULL)||(F2==NULL)||(F3==NULL))
{
printf("Erreur d ouverture d un ou des fichiers texte de resultats\n");
getch();
exit(1);
}
/*ecriture des 18 lignes de l'entete */
fprintf(F," calcul des profils d'aile par la methode des panneaux\n");
fprintf(F,"méthode avec tourbillons par formulation de ligne de courant\n\n");
fprintf(F,"BOLOVIS Angela FOISSAC Émilie LEROUYER Benjamin MEUNIER Arnaud\n\n");
if (naca < 100) fprintf (F,"Les caracteristiques du profil NACA symetrique 00%d sont les suivantes :\n",naca);
else fprintf (F,"Les caracteristiques du profil NACA cambre %d sont les suivantes :\n",naca);
fprintf (F," * Longueur de reference : c = 1 pour l etude\n");
fprintf (F," * Position de la cambrure maximale : d = %lf\n",*d);
fprintf (F," * Cambrure maximale : f = %lf\n", *f);
fprintf (F," * Rayon du bord d'attaque : r = %lf\n", *r);
fprintf (F,"\n\n L angle d incidence est de %lf degres.\n", alpha);
fprintf (F," Pour %d panneaux\n", N);
fprintf(F,"le coefficient de portance du profile est : %lf\n\n\n",Cl[N]);
fprintf(F,"X PtCont | Y PtContr | CoefPress | CoefPorta | X Noeud | Y Noeud\n");
/* ecriture des N enregistrements portant les valeurs caracteristique des panneaux i */
for(i=0;i<N;i++)
{
fprintf(F,"%lf | %+lf | %+lf | %+lf | %lf | %+lf\n",X[i],Y[i],Cp[i], Cl[i], Xi[i], Yi[i]);
fprintf(F1,"%lf %+lf %+lf %+lf %lf %+lf\n",X[i],Y[i],Cp[i], Cl[i], Xi[i], Yi[i]);
fprintf(F2,"%lf %+lf\n",X[i],Cp[i]);
fprintf(F3,"%lf %+lf %+lf\n",X[i],Y[i],Cp[i]);
}
fprintf(F1,"%lf %+lf %+lf %+lf %lf %+lf\n",X[N],Y[N],Cp[0], Cl[N] ,Xi[N], Yi[N]);
/* fermeture des fichiers resultat.txt et resultataero.txt */
fclose(F);
fclose(F1);
fclose(F2);
fclose(F3);
}/*fin de la fonction ecrire fichier */
/****************************************************************************/
/* ecriture du fichier apres un retraitement du programme pour le meme NACA */
/****************************************************************************/
void ecrirefichier1(double *X,double *Y,double *Xi,double *Yi,double *Cl,double *Cp, double alpha, int choix)
{
int i;
FILE *F=fopen("resultataero.txt","a+");
FILE *F1=fopen("resultat.txt","a+");
FILE *F2=fopen("resultat_X_Cp.txt","a+");
FILE *F3=fopen("resultat_X_Y_Cp.txt","a+");
/* si ouverture du fichier impossible*/
if ((F==NULL)||(F1==NULL)||(F2==NULL)||(F3==NULL))
{
printf("Erreur d ouverture d un ou des fichiers texte de resultats\n");
getch();
exit(1);
}
fprintf(F,"\n\n--- nouveau traitement avec le nombre de panneau porte a %d ---\n", N);
if (choix==3) fprintf(F,"--- et un nouvel angle d incidence alpha = %lf ---\n", alpha);
fprintf(F,"\n");
fprintf(F2,"\n");
fprintf(F3,"\n");
fprintf(F,"X PtContr | Y PtCont | CoefPress | CoefPorta | X Noeud | Y Noeud\n");
/* ecriture des N enregistrements portant les valeurs caracteristique des panneaux i */
for(i=0;i<N;i++)
{
fprintf(F,"%lf | %+lf | %+lf | %+lf | %lf | %+lf\n",X[i],Y[i],Cp[i], Cl[i], Xi[i], Yi[i]);
fprintf(F1,"%lf %+lf %+lf %+lf %lf %+lf\n",X[i],Y[i],Cp[i], Cl[i], Xi[i], Yi[i]);
fprintf(F2,"%lf %+lf\n",X[i],Cp[i]);
fprintf(F3,"%lf %+lf %+lf\n",X[i],Y[i],Cp[i]);
}
fprintf(F1,"%lf %+lf %+lf %+lf %lf %+lf\n", X[N],Y[N], Cp[0], Cl[N], Xi[N], Yi[N]);
/* fermeture des fichiers resultat.txt et resultataero.txt */
fclose(F);
fclose(F1);
fclose(F2);
fclose(F3);
}/*fin de la fonction ecrire fichier deuxieme passage */
/**********************************************************************************************/
/* ecriture du fichier du coefficient de portance Cl en fonction de l angle d incidence alpha */
/**********************************************************************************************/
void ecrirefichier2(double *Cl, double alpha, double alpha2, int naca, double *d, double *f, double *r, int *spec1)
{
FILE *F;
FILE *F1;
/* ouverture pour la premiere fois du fichier: on efface le precedent et on imprime l'entête*/
if (*spec1==0)
{
*spec1 = 1;
F=fopen("resultat_alpha_Cl_aero.txt","w");
F1=fopen("resultat_alpha_Cl.txt","w");
/* si ouverture du fichier impossible*/
if ((F==NULL)||(F1==NULL))
{
printf("Erreur d ouverture du fichier resultat_alpha_Cl\n");
getch();
exit(1);
}
fprintf(F,"Fichier resultat de calcul du coefficient de portance Cl en fonction de l angle d incidence alpha\n");
fprintf(F,"méthode avec tourbillons par formulation de ligne de courant\n\n");
fprintf(F,"BOLOVIS Angela FOISSAC Emilie LEROUYER Benjamin MEUNIER Arnaud\n\n");
if (naca < 100) fprintf (F,"Les caracteristiques du profil NACA symetrique 00%d sont les suivantes :\n",naca);
else fprintf (F,"Les caracteristiques du profil NACA cambre %d sont les suivantes :\n",naca);
fprintf (F," * Longueur de reference : c = 1 pour l etude\n");
fprintf (F," * Position de la cambrure maximale : d = %lf\n",*d);
fprintf (F," * Cambrure maximale : f = %lf\n", *f);
fprintf (F," * Rayon du bord d'attaque : r = %lf\n", *r);
fprintf (F,"\n\n l angle d incidence varie de %lf a %lf degres.\n", alpha, alpha2);
fprintf (F," Pour %d panneaux\n\n", N);
fprintf(F,"alpha | Cl du profil\n");
}
else
{
F=fopen("resultat_alpha_Cl_aero.txt","a+");
F1=fopen("resultat_alpha_Cl.txt","a+");
/* si ouverture du fichier impossible*/
if ((F==NULL)||(F1==NULL))
{
printf("Erreur d ouverture du fichier resultat_alpha_Cl\n");
getch();
exit(1);
}
}
/* ecriture de l enregistrement */
fprintf(F,"%lf %+lf\n",alpha,Cl[N]);
fprintf(F1,"%lf %+lf\n",alpha,Cl[N]);
/* fermeture des fichiers resultat.txt et resultataero.txt */
fclose(F);
fclose(F1);
}/* fin de la fonction ecrire fichier Cl fonction de alpha */
/**************************************************************************/
/* fonction principale main */
/**************************************************************************/
void main ()
{
int naca, choix=0;
double d, f, r, t, v_inf, alpha, alpha2=0;
double **A, *Xi, *Yi, *X, *Y;
double *l, *b, *Cp, *Cl;
char sym;
/* cet indicateur sert à savoir si l'on a déjà ouvert notre fichier(=1) ou non */
int spec1 = 0;
/* Saisie et calcul des caracteristiques du profil */
menu(&naca, &v_inf, &alpha);
caracteristiques (naca, &sym, &d, &f, &r, &t);
while (choix!= 1)
{
/* allocation des vecteurs contenant la position des noeuds */
Xi = alloc_vec (Xi,N+1);
Yi = alloc_vec (Yi,N+1);
/* discretisation du profil */
discretisationNACA (sym, f, d, t, Xi, Yi);
/* Allocation des matrices et vecteurs remplis par la fonction initmatrice*/
A = alloc_mat (A,N+1,N+1);
l = alloc_vec (l,N);
/* on alloue une position supplementaire car le premier element et le dernier
sont les memes */
X = alloc_vec (X,N+1);
Y = alloc_vec (Y,N+1);
/* remplissage de la matrice A des coefficients */
initmatrice (A, Xi, Yi, X, Y, l);
b = alloc_vec (b, N+1);
/* remplissage du vecteur b des vitesses */
init_b (b, &v_inf, alpha, X, Y);
/* appel a la fonction gauss. On envoie le vecteur b, et l on recupere
en sortie le vecteur x dans le vecteur b */
if (gauss_pivmax(A,b)==0)
{
printf("erreur dans la resolution du systeme");
getch();
exit(1);
}
/* calcul des coefficients de pression et de portance */
Cp = alloc_vec (Cp, N);
C_pression (b, &v_inf, Cp);
/* la derniere position du vecteur portance contient la portance totale
sur l'aile */
Cl = alloc_vec (Cl, N+1);
C_portance (b, &v_inf, Cl,l);
/*ecriture du fichier resultat en sortie nous permettant de tracer
des graphiques sous MAPLE ou excel */
if (choix == 0) ecrirefichier(X, Y, Xi, Yi, Cl, Cp, naca, &d, &f, &r, alpha);
if ((choix == 2)||(choix == 3)) ecrirefichier1(X, Y, Xi, Yi, Cl, Cp, alpha, choix);
if (choix == 4) ecrirefichier2(Cl, alpha, alpha2, naca, &d, &f, &r, &spec1);
/* liberation de la memoire vive pour reprendre ensuite des matrices et des vecteurs de tailles differentes
si N est different */
free(A);
free(l);
free(Xi);
free(Yi);
free(Cp);
free(Cl);
/* on recommence le programme pour une autre incidence ou pour une courbe des Cl en fonction des alphas */
if ((choix==4)&&(alpha<(alpha2))) alpha++;
else choix = menu2(naca, &alpha, &alpha2);
} /* fin du while de reprise des traitements */
printf ("\nle programme s est deroule correctement.\n Vous pouvez voir les resultats dans le fichier resultataero.txt\n");
printf ("\nLe fichier resultat.txt charge sous excel permet de visualiser les graphiques\n");
getch();
}/* fin de la fonction main */