基于中心差分有限离散化和 Newton Raphson 算法求解NACA 翼型二维不可压缩和可压缩流动附matlab代码

✅作者简介：热爱科研的Matlab仿真开发者，修心和技术同步精进，matlab项目合作可私信。

🍎个人主页：Matlab科研工作室

🍊个人信条：格物致知。

更多Matlab仿真内容点击👇

智能优化算法       神经网络预测       雷达通信       无线传感器        电力系统

信号处理              图像处理               路径规划       元胞自动机        无人机 

⛄ 内容介绍

该程序基于使用小微扰理论的方法求解组合欧拉方程，使用中心差分有限离散化和牛顿拉斐尔森算法。

⛄ 部分代码

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Programme principal

clear all;    % Effacer les variables

clc;          % Effacer l'ecran

format long   % Affichage en double précision 

%% Données du code

prompt={'Maximum iterations', 'Calculation accuracy','Domain length (m)', ...

    'Domain Width (m)','X Step', 'Y Step','Start Airfoil position (m)',...

    'Airfoil chorde length (m)','Initial estimatation','Mach number','Airfoil thickness'};

title='Airfoils Calculator';

lines=1;

def={'100','1e-06','4.0','2.5','0.05','0.05','1.5','1.0','0.5','0.1','0.006'};

answer=inputdlg(prompt,title,lines,def);

assignin('base','maxIter',str2double(answer{1})); % max des iterations

assignin('base','tolX',str2double(answer{2}));% erreur de calcul

assignin('base','Lx',str2double(answer{3}));% longeur du domaine 

assignin('base','Ly',str2double(answer{4})); % largeur du domaine

assignin('base','hx',str2double(answer{5})); % pas suivant x

assignin('base','hy',str2double(answer{6})); % pas suivant y

assignin('base','Dp',str2double(answer{7})); % position de début de profil

assignin('base','C',str2double(answer{8})); % la corde du profil

assignin('base','est_x',str2double(answer{9})); % Estimation initiale

assignin('base','M',str2double(answer{10})); % Mach number

assignin('base','t',str2double(answer{11})); % Airfoil thickness

%%

Nx=ceil(Lx/hx);  % nombre des éléments suivant ox

Ny=ceil(Ly/hy);  % nombre des éléments suivant oy

N=Nx*Ny; % nombre des équations à résoudre

Fp=Dp+C; % calcul la position de la fin du profil

Dp=ceil(Dp/hx); % position de début du maillage

Fp=ceil(Fp/hx);  % position de la fin du maillage

%% Pour tracer graphiquement les résultats

debprof=Dp:Fp;

debx=2:Nx-1;

deby=2:Ny-1;

%% Génération du maillage 

x_grid=linspace(0,Lx,Nx+1);

y_grid=linspace(0,Ly,Ny+1);

[X1,Y1]=meshgrid(x_grid,y_grid);

[X2,Y2]=meshgrid(y_grid,x_grid);

%% % Conditions initiales

for i=1:N

X0(i)=est_x;   % Estimation Initiale

end

X0=X0';        % X0 sous forme colonne 

%% Méthode de Newton Raphson

x=X0;

Xoled=X0;

disp('_____________________________Iterations______________________________')

for i=1:maxIter

    F=EQUATIONS(x,Dp,Fp,Nx,Ny,hx,hy,N,M,t);

    jac=JACFUNC(F,x,Dp,Fp,Nx,Ny,hx,hy,N,M,t);

    x=x-(inv(jac))*F;

    err(:,i)=abs(x-Xoled);

    Xoled=x;

    if (err(:,i)< tolX)     

      break;

    end 

    disp(['Iter   ',num2str(i),'   Err= ',num2str(max(err(:,i)))]);

Err(i)=max(err(:,i));

iter(i)=i;

end

disp('_______________________________Slutions______________________________')

%% Affichage les solutions

for i=1:N

disp(['x(',num2str(i),')=  ',num2str(x(i))])

end

%% Transformation vecteur x en Matrice Phi

cont=0;

for j=1:Ny

    for i=1:Nx

phi(i,j)=x(i+Nx*cont);

    end 

    cont=cont+1;

end

%% Calcul - vitesse U

for i=2:Nx-1

for j=1:Ny

U(i,j)=1+((phi(i+1,j)-phi(i-1,j)))/(2*hx);

end

end

%% Calcul - vitesse V 

for i=2:Nx-1

for j=2:Ny-1

V(i,j)=((phi(i,j+1)-phi(i,j-1)))/(2*hy);

end

end

%% Calcul - Coefficient de pression (Cp)

for i=2:Nx-1

    for j=1:Ny

   CP(i,j)=-2*((phi(i+1,j)-phi(i-1,j)))/(2*hx);

    end

end

%% calcul Cn et Ct

Cn=0;

for i=Dp:Fp

    for j=1:2

Cn=((CP(i,j)+CP(i+1,j))/2)*((i+1)*hx-i*hx)+Cn;

end

end

%

Ct=0;

for i=Dp:Fp

    for j=1:2

Ct=((CP(i,j)+CP(i+1,j))/2)*((j+1)*hy-j*hy)+Ct;

    end

end

%% Calcul Cl et Cdp

Alpha=0.0;

Cl=Cn*sin(Alpha)+Ct*cos(Alpha);

Cdp=Cn*sin(Alpha)-Ct*cos(Alpha);

disp('____________________________________________________________________')

%% Affichage graphique

figure

subplot(2,2,1)

box on

grid on

plot(iter,Err,'ro-')

xlabel('Iter')

ylabel('Err')

%title('Convergence Courbe')

%% Afficher les figures 

subplot(2,2,2)

hold on 

grid on 

box on 

plot(phi(:,1))

plot(phi(:,2))

plot(phi(:,3))

plot(phi(:,4))

plot(phi(:,5))

plot(phi(:,6))

plot(phi(:,7))

plot(phi(:,Ny))

legend('y=1','y=2','y=3','y=4','y=5','y=6','y=7','y=Ny')

%title('Fonction de courant Phi')

xlabel('X (m)')

ylabel('Phi')

hold off

%% Tracer la vitesse U 

subplot(2,2,3)

hold on 

grid on 

box on 

plot(U(debx,1))

plot(U(debx,2))

plot(U(debx,3))

plot(U(debx,4))

plot(U(debx,5))

plot(U(debx,6))

plot(U(debx,7))

plot(U(debx,Ny-1))

legend('y=1','y=2','y=3','y=4','y=5','y=6','y=7','y=Ny-1')

%title('Vitesse U/Uinf')

xlabel('X (m)')

ylabel('U/Uinf')

hold off

%% Tracer la vitesse V 

subplot(2,2,4)

hold on 

grid on 

box on 

debn=15:25;

%plot(V(deb,1))

plot(V(:,2))

plot(V(:,3))

plot(V(:,4))

plot(V(:,5))

plot(V(:,6))

plot(V(:,7))

plot(V(:,8))

legend('y=2','y=3','Y=4','y=5','y=6','y=7','Y=8')

%title('Vitesse V/Vinf')

xlabel('X (m)')

ylabel('V/Vinf')

hold off

%% Tracer CP

figure

subplot(2,2,1)

hold on 

grid on 

box on 

%plot(V(deb,1))

plot(-CP(debprof,2))

plot(-CP(debprof,3))

plot(-CP(debprof,4))

plot(-CP(debprof,5))

plot(-CP(debprof,6))

plot(-CP(debprof,7))

plot(-CP(debprof,8))

legend('y=2','y=3','Y=4','y=5','y=6','y=7','Y=8')

%title(' Coefficient de pression')

xlabel('X (m)')

ylabel('-CP')

hold off

%% Afficher les contours 

subplot(2,2,2)

[C,h]=contourf(phi,25,'ShowText','on')

xlabel('Contour Phi')

subplot(2,2,3)

contourf(U(debx,deby),25)

xlabel('Contour U')

subplot(2,2,4)

contourf(V(debx,deby),25)

xlabel('Contour V')

figure

subplot(2,1,1)

contourf(CP,25)

xlabel('Contour CP')

subplot(2,1,2)

hold on

plot(X1,Y1,'k-',Y2,X2,'k-')

axis([0 Lx 0 Ly])

%title('Mesh Grid')

box on

hold off

%_____________________________________ end _________________________________

⛄ 运行结果

⛄ 参考文献

⛳️ 完整代码

❤️部分理论引用网络文献，若有侵权联系博主删除

❤️ 关注我领取海量matlab电子书和数学建模资料