Implémentation TDMA en 2D en C (Mécanique des fluides)
Répondre à la discussion
Affichage des résultats 1 à 5 sur 5

Implémentation TDMA en 2D en C (Mécanique des fluides)



  1. 05/03/2017, 23h33 #1
    Tash

    Implémentation TDMA en 2D en C (Mécanique des fluides)


    ------

    Bonjour,

    Je souhaiterais implémenter la méthode TDMA(algo de Thomas) en 2D, mais ya un truc qui cloche dans mon programme :

    Code:
    #include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<string.h>
#define Jmax 11
#define Imax 12
int main()
{
	int i,it,l,k,j;
	//****les constantes
	double nu=1e-3;
	double	dy=0.02;
	//le fichier de sortie(pour gnuplot)
	FILE *f;
	//les tableaux
	double a[Imax],b[Imax],c[Imax],d[Imax];
	double Ue[Imax];
	double x[Imax];
	double dx[Imax];
	//***************
	double P[Imax];
	double Q[Imax];
		// on crée les tableaux à deux dimensions U et V
		double U[Imax][Jmax];
		double V[Imax][Jmax];
//*****************************************************
	for(i=1;i<Imax;i++){
		x[i]=i+0.5;
	}
	for(i=1;i<=Imax;i++){
		dx[i]=x[i+1]-x[i]/2;
	}
	//
	for(i=1;i<Imax;i++)
	{
		Ue[i]=sin(x[i]);
	}
	//
	//Conditions aux limites***********************************************************
	//bord gauche
	for(i=1;i<=Imax;i++)
	{
		U[i][1]=0.000;
		V[i][1]=0.000;
	}
	//bord droite
	for(i=1;i<=Imax;i++)
	{
		U[i][Jmax]=Ue[i];
		V[i][Jmax]=0.0;
	}
	//bord extérieur(en haut du maillage)
	for(j=1;j<Jmax;j++)
	{
		U[1][j]=0.000; 
		V[1][j]=0.000;
	}
	//*******************************************************************************************
	//**************************Application : Méthode TDMA***************************************
	//*******************************************************************************************
	//*******sur les lignes Nord-Sud ************************************************************
	//*******************************************************************************************
	P[1]=0.0;
	for(i=2;i<Imax;i++)
	{	
		Q[1]=U[i][1];   		//tous nuls(Conditions aux limites gauche)
		for(j=2;j<Jmax;j++)
		{
			//calcul des coef. ligne par ligne
			a[j]=(V[i][j]/2*dy)-(nu/(dy*dy));
			b[j]=-(V[i][j]/2*dy)-(nu/(dy*dy));
			d[j]=(U[i][j]/dx[i])+(2*nu/(dy*dy));
			c[j]=(U[i][j]*U[i][j])/dx[i]+(Ue[i]*(Ue[i+1]-Ue[i])/dx[i]);
			//**********************************************************
			//on Calcule P(.,j) et Q(.,j) pour chaque ligne de j=2:Jmax
			//**********************************************************
			P[j]=-a[j]/(d[j]+b[j]*P[j-1]);
			Q[j]=(c[j]-b[j]*Q[j-1])/(d[j]+b[j]*P[j-1]);
			//**********************************************************
			//puis on calcule  U(.,j) pour chaque ligne de j=Jmax-1 à 1
			//**********************************************************
			do{
				U[i][j]=P[j]*U[i][j-1]+ Q[j];
				j=j-1;
			}while(j!=1);
		}
		//**************************************************************
		//une fois U(.,j) obtenu, on calcule V(.j) j=1:Jmax-1
		//**************************************************************
		for(j=1;j<Jmax;j++)
		{	V[i][j]=V[i][j-1]-(dy/dx[i])*(U[i][j-1]-U[i][j])+(dy/2.*dx[i])*(U[i][j-1]-U[i][j]);}
	}//*********************fin Nord-Sud*********************************************************
	//*******************************************************************************************
	//********************sur les colonnes Est-Ouest*********************************************
		P[1]=0.0;
	for(j=2;i<Jmax;j++)
	{	
		Q[1]=U[1][j];   		//tous nuls(Conditions aux limites gauche)
		for(i=2;i<Imax;i++)
		{
			//calcul des coef. colonne par colonne
			a[i]=(V[i][j]/2*dy)-(nu/(dy*dy));
			b[i]=-(V[i][j]/2*dy)-(nu/(dy*dy));
			d[i]=(U[i][j]/dx[j])+(2*nu/(dy*dy));
			c[i]=(U[i][j]*U[i][j])/dx[j]+(Ue[j]*(Ue[j+1]-Ue[j])/dx[j]);
			//**********************************************************
			//on Calcule P(i,.) et Q(i,.) pour chaque colonne de i=2:Imax-1
			//**********************************************************
			P[i]=-a[i]/(d[i]+b[i]*P[i-1]);
			Q[i]=(c[i]-b[i]*Q[i-1])/(d[i]+b[i]*P[i-1]);
			//**********************************************************
			//puis on calcule  U(i,.) pour chaque colonne de i=Imax-1:2
			//**********************************************************
			do{
				U[i][j]=P[i]*U[i-1][j]+ Q[i];
				printf("U[%d][%d]\n",i,j,U[i][j]);
				i=i-1;
			}while(i!=2);
		}
		//**************************************************************
		//une fois U(i,.) obtenu, on calcule V(i,.) i=1:Imax-1
		//**************************************************************
		for(i=1;i<Imax;i++)
		{	V[i][j]=V[i-1][j]-(dy/dx[j])*(U[i-1][j]-U[i][j])+(dy/2.*dx[j])*(U[i-1][j]-U[i][j]);
			printf("V[%d][%d]\n",i,j,V[i][j]);
		}
	}//*********************fin Est-Ouest********************************************************
	//************************************************************************************************************
	//**************************FIN TDMA-2D(U,V)******************************************************************
	//************************************************************************************************************
	//******************écriture des données pour gnuplot****************************************
	//**********************je veux simuler le profil de vitesse(i.e U/Ue en fonction x/c=0.5;1;1.5;2;2.5,...)*****
			f=fopen("meca-result.dat","w");
			for(i=1;i<=Imax;i++)
			{
				for(j=1;j<=Jmax;j++)
				{
					fprintf(f,"%lf %lf \n",U[i][j]/Ue[i],x[i]);
				}
			}
			fclose(f);
	return 1;
}

    -----
    Dernière modification par JPL ; 06/03/2017 à 00h02. Motif: ajout de la balise Code pour garder l'indentation
  2. 05/03/2017, 23h50 #2
    Chanur

    Re : Implémentation TDMA en 2D en C (Mécanique des fluides)

    Bonjour,

    Si vraiment tu veux de l'aide, ce serait aussi bien de ne pas faire tout ton possible pour compliquer les choses ...
    => utilise de balises CODE : là, sans les indentations c'est carrément illisible.
    => "y a un truc qui cloche" : c'est vraiment l'indication la plus précise que tu puisses donner ? tu ne sais même pas si ton programme se compile ? s'il plante à l'exécution ? s'il fournit des résultats faux ? ou mal formatés ? s'il provoque des fuites radio-actives ? Parce que moi (et tous les autres lecteurs) j'en suis réduit à essayer de le deviner ...
    => tu sais qu'en langage C, on a le droit de donner aux variables un nom qui a un rapport avec leur signification ? Parce que a,b,c,d,i,j,k,l,x etc. c'est carrément cryptique ...
    Ce qui se conçoit bien s'énonce clairement ; et les mots pour le dire arrivent aisément.
  3. 06/03/2017, 00h05 #3
    Chanur

    Re : Implémentation TDMA en 2D en C (Mécanique des fluides)

    Bon, mon coup de colère passé, et comme mon compilateur n'a pas d'états d'âmes, il m'a suffit d'essayer de faire tourner ton programme sous debugger pour y trouver une boucle sans fin :
    Code:
    for(j=2;j<Jmax;j++) // j'incrémente j (un pas en avant)
   {
    //calcul des coef. ligne par ligne
    a[j]=(V[i][j]/2*dy)-(nu/(dy*dy));
    b[j]=-(V[i][j]/2*dy)-(nu/(dy*dy));
    d[j]=(U[i][j]/dx[i])+(2*nu/(dy*dy));
    c[j]=(U[i][j]*U[i][j])/dx[i]+(Ue[i]*(Ue[i+1]-Ue[i])/dx[i]);
    //****************************** ****************************
    //on Calcule P(.,j) et Q(.,j) pour chaque ligne de j=2:Jmax
    //****************************** ****************************
    P[j]=-a[j]/(d[j]+b[j]*P[j-1]);
    Q[j]=(c[j]-b[j]*Q[j-1])/(d[j]+b[j]*P[j-1]);
    //****************************** ****************************
    //puis on calcule U(.,j) pour chaque ligne de j=Jmax-1 à 1
    //****************************** ****************************
    do
       {
         U[i][j]=P[j]*U[i][j-1]+ Q[j];
         j=j-1; // je fais redescendre j à 1 (un pas en arrière)
       }while(j!=1);
   } // tant que j ne vaut pas 11 : On n'est pas arrivé ...
    Tu as un debugger ?
    Ce qui se conçoit bien s'énonce clairement ; et les mots pour le dire arrivent aisément.
  4. 06/03/2017, 11h18 #4
    Tash

    Re : Implémentation TDMA en 2D en C (Mécanique des fluides)

    Bonjour,
    mes excuses pour la présentation, je suis nouveau et j'ai pas pris le temps de jeter un coup d'oeil sur le réglement(j'ai eu tort!).

    en fait, c'est un système d'équations en 2D avec deux composantes(U,V),,
    et t'as pas tort, j'ai du mal copier-coller(la boucle for sur j dois être fermée avant l'entrée du boucle dou--while,,,

    i pour parcourir l'axe horizontal, j pour l'axe vertical

    le système à résoudre :
    Code:
    $$\left\{ 
\begin{array}{cc}\label{m6}
\hspace{0.5cm}U_i^{j+1}=P_iU_{i+1}^{j+1}+Q_i\quad i=2:Imax-1,\dots,j=2:Jmax-1\vspace{0.3cm}\\
V_i^{j+1}=V_{i-1}^{j+1}-e_i(U_{i-1}^{j+1}-U_i^{j+1})+d_j\quad\quad i=1:Imax,\dots,j=1:Jmax
\end{array}\right.$$
Avec 
$$\left\{ 
\begin{array}{cc}\label{m6}
P_i=-A_i\left(D_i+B_iP_{i-1}\right)^{-1}\vspace{0.3cm}\nonumber\\
Q_i=\left(C_i-B_iQ_{i-1}\right)\left(D_i+B_iP_{i-1}\right)^{-1}\nonumber\vspace{0.3cm}\\
B_i=-\frac{V_i^j}{2\Delta y}-\frac{\nu}{(\Delta y)^2}\nonumber\vspace{0.3cm}\\
D_i=\frac{U_i^j}{\Delta x}+\frac{2\nu}{(\Delta y)^2}\nonumber\vspace{0.3cm}\\
A_i=\frac{V_i^j}{2\Delta y}-\frac{\nu}{(\Delta y)^2}\nonumber\vspace{0.3cm}\\
C_i=\frac{(U_i^j)^2}{\Delta x}+U_e^j\frac{(U_e^{j+1}-U_e^j)}{\Delta x}\nonumber\vspace{0.3cm}\\
d_j=\frac{\Delta y}{2\Delta x}(U_{i-1}^j+U_i^j)\nonumber\vspace{0.3cm}\\
e_i=\frac{\Delta y}{\Delta x}\nonumber\vspace{0.3cm}\\
P_1=Q_1=0\quad\quad P_{Imax}=0\quad Q_{Imax}=U_e\quad U_e<0.99) \nonumber\vspace{0.3cm}\\
V_1^{j+1}=0,\quad V_{Imax}^{j+1}=0\nonumber\vspace{0.3cm}\\
\end{array}\right.$$
Pour l'instant,j'ai choisi $x=0.5,1,1.5,\dots$, et $dx=x_{i+1}-xi/2$, $Ue=sin(x)$ $dy=0.02$ $U_1=0$ et $U_N=Ue$, $nu=1e-3$
    le principe de l'algo :
    Code:
    -Calculer (P_2,Q_2),(P_3,Q_3),...(P_{N-1},Q_{N-1})
_puis U_{N-1},U_{N-2},...,U_2
_puis V_1,...,V_{N-1}(à partir de U)
    Code:
    *********************************************************************************************
rectification:(sur la ligne Nord-Sud)
//*******sur les lignes Nord-Sud ************************************************************
	//*******************************************************************************************
	//*******sur les lignes Nord-Sud ************************************************************
	//*******************************************************************************************

	P[1]=0.0;
	for(i=2;i<Imax;i++)
	{	
		Q[1]=U[i][1];   		//tous nuls(Conditions aux limites gauche)
		for(j=2;j<Jmax;j++)
		{
			//calcul des coef. ligne par ligne
			a[j]=(V[i][j]/2*dy)-(nu/(dy*dy));
			b[j]=-(V[i][j]/2*dy)-(nu/(dy*dy));
			d[j]=(U[i][j]/dx[i])+(2*nu/(dy*dy));
			c[j]=(U[i][j]*U[i][j])/dx[i]+(Ue[i]*(Ue[i+1]-Ue[i])/dx[i]);
			//
			//**********************************************************
			//on Calcule P(.,j) et Q(.,j) pour chaque ligne de j=2:Jmax
			//**********************************************************
			P[j]=-a[j]/(d[j]+b[j]*P[j-1]);
			Q[j]=(c[j]-b[j]*Q[j-1])/(d[j]+b[j]*P[j-1]);
			//**********************************************************
			//puis on calcule  U(.,j) pour chaque ligne de j=Jmax-1 à 1
			//**********************************************************
			k=Jmax-1 ;
			do{
				U[i][k]=P[j]*U[i][k+1]+ Q[j];
				k=k-1;
			}while(k!=2);
		//**************************************************************
		//une fois U(.,j) obtenu, on calcule V(.j) j=1:Jmax-1
		//**************************************************************
		for(l=1;l<Jmax;l++)
		{	V[i][l]=V[i][l-1]-(dy/dx[i])*(U[i][l-1]-U[i][l])+(dy/2.*dx[i])*(U[i][l-1]-U[i][l]);}
	}
}//*********************fin Nord-Sud************
	//*******************************************************************************************
  5. Aujourd'hui
    A voir en vidéo sur Futura
  6. 06/03/2017, 11h22 #5
    Tash

    Re : Implémentation TDMA en 2D en C (Mécanique des fluides)

    ci-joint le sujet en pdf ...
    Images attachées Images attachées
« Makefile sur Windows 7 avec Gcc | template c++ »

Discussions similaires

  1. Mécanique des fluides
    Par mango1 dans le forum Physique
    Réponses: 8
    Dernier message: 29/06/2015, 14h12
  2. Tdma
    Par marwa8888 dans le forum Électronique
    Réponses: 1
    Dernier message: 25/12/2014, 12h02
  3. Mecanique des fluides
    Par Sawyerbond dans le forum Physique
    Réponses: 3
    Dernier message: 21/12/2013, 06h50
  4. Mécanique des fluides
    Par bboop8 dans le forum Physique
    Réponses: 3
    Dernier message: 31/10/2008, 14h50
  5. Mecanique des fluides et statique des fluides
    Par invitef1754d56 dans le forum Physique
    Réponses: 0
    Dernier message: 09/12/2007, 12h25