该系列MATLAB代码实现了一种变稳方法（VortexLatticeMethod,VLM），用于计算机翼的气动特性。主要函数包括`VLM.m`、`mesh.m`、`VLMalphaList.m`等，涉及网格生成、涡旋分布计算、气动力计算等步骤。代码首先定义了不同角度攻角（alpha_deg）和网格数（nCells）的组合，然后对每个组合进行VLM计算并绘制结果，展示了不同攻角和网格数量对机翼升力和阻力的影响。
main.m
clearvars; clc; % main.m
f_max=0.08; % f_max=0.0;
alpha_deg_list = [-10;-5;0;5;10]; nCellsList = [8;16;32;64;128];
objVLMalphaList = VLMalphaList(f_max,alpha_deg_list,nCellsList);
solve_plot(objVLMalphaList);
mesh.m
% mesh.m
classdef mesh
properties
nCells
f_max
end
methods
function obj = mesh(nCells,f_max)
obj.nCells = nCells;
obj.f_max = f_max;
end
function y = get_y(obj,x)
c = 1;
y = obj.f_max*(1-(2*x/c-1)^2);
end
function xy_mat = mesh_points(obj)
c = 1;
delta_x = c / obj.nCells;
x = zeros(obj.nCells+1,1);
y = zeros(size(x,1),1);
for i = 1 : size(x,1)
x(i,1) = (i-1)*delta_x;
y(i,1) = get_y(obj,x(i,1));
end
xy_mat = [x,y];
end
function xy_mat = mesh_slope(obj)
xy_p = mesh_points(obj);
x_p = xy_p(:,1);
y_p = xy_p(:,2);
y_k = zeros(obj.nCells,1);
x_k = zeros(size(y_k,1),1);
for i = 1 : size(x_k,1)
y_k(i,1) = (y_p(i+1,1)-y_p(i,1))/(x_p(i+1,1)-x_p(i,1));
x_k(i,1) = (x_p(i+1,1)+x_p(i,1))/2;
end
xy_mat = [x_k,y_k];
end
function xy_mat = mesh_linear_points(obj,fraction)
xy_p = mesh_points(obj);
xy_k = mesh_slope(obj);
y_k = xy_k(:,2);
x_p = xy_p(:,1);
y_p = xy_p(:,2);
x_line = zeros(obj.nCells,1);
y_line = zeros(size(x_line,1),1);
for i = 1 : size(y_line,1)
delta_x = x_p(i+1,1) - x_p(i,1);
delta_x = delta_x*fraction;
x_line(i,1) = x_p(i,1) + delta_x;
y_line(i,1) = y_p(i,1) + delta_x*y_k(i,1);
end
xy_mat = [x_line,y_line];
end
function xy_mat = mesh_vortex_points(obj)
fraction = 0.25;
xy_mat = mesh_linear_points(obj,fraction);
end
function xy_mat = mesh_middle_points(obj)
fraction = 0.5;
xy_mat = mesh_linear_points(obj,fraction);
end
function xy_mat = mesh_control_points(obj)
fraction = 0.75;
xy_mat = mesh_linear_points(obj,fraction);
end
function xyr_mat = mesh_vortex_xyr(obj,xp,yp)
x_mat = zeros(obj.nCells,1);
y_mat = zeros(size(x_mat,1),1);
r_mat = zeros(size(y_mat,1),1);
xy_v = mesh_vortex_points(obj);
x_v = xy_v(:,1);
y_v = xy_v(:,2);
for i = 1 : size(r_mat,1)
xx = xp - x_v(i,1);
yy = yp - y_v(i,1);
r_mat(i,1) = sqrt(xx^2+yy^2);
x_mat(i,1) = yy;
y_mat(i,1) = - xx;
end
xyr_mat = [x_mat,y_mat,r_mat];
end
function deltax = get_dx(obj)
c = 1;
deltax = c / obj.nCells;
end
function deltac = get_dc(obj,ii)
xy_p = mesh_points(obj);
x_p = xy_p(:,1);
y_p = xy_p(:,2);
dx = x_p(ii+1,1) - x_p(ii,1);
dy = y_p(ii+1,1) - y_p(ii,1);
deltac = sqrt(dx^2+dy^2);
end
end
end
VLM.m
% VLM.m
classdef VLM
properties
f_max
alpha
nCells
end
methods
function obj = VLM(f_max,alpha_deg,nCells)
obj.f_max = f_max;
obj.alpha = alpha_deg / 180 * pi;
obj.nCells = nCells;
end
function plot_mat(~,mat)
x = mat(:,1);
y = mat(:,2);
plot(x,y);
end
function scatter_mat(~,mat)
x = mat(:,1);
y = mat(:,2);
scatter(x,y);
end
function xy_bb = get_bb(obj)
U = 1;
y_b = zeros(obj.nCells,1);
Uy = U * sin(obj.alpha);
Ux = U * cos(obj.alpha);
objMesh = mesh(obj.nCells,obj.f_max);
xy_k = mesh_slope(objMesh);
xy_c = mesh_control_points(objMesh);
x_c = xy_c(:,1);
y_k = xy_k(:,2);
for i = 1 : size(y_b,1)
y_b(i,1) = Uy - y_k(i,1)*Ux;
end
xy_bb = [x_c,y_b];
end
function A_mat = get_AA(obj)
A_mat = zeros(obj.nCells,obj.nCells);
objMesh = mesh(obj.nCells,obj.f_max);
xy_c = mesh_control_points(objMesh);
xy_k = mesh_slope(objMesh);
y_k = xy_k(:,2);
x_c = xy_c(:,1);
y_c = xy_c(:,2);
for ii = 1 : size(A_mat,1)
xp = x_c(ii,1);
yp = y_c(ii,1);
xyr_v = mesh_vortex_xyr(objMesh,xp,yp);
x_v = xyr_v(:,1);
y_v = xyr_v(:,2);
r_v = xyr_v(:,3);
for jj = 1 : size(A_mat,2)
A_mat(ii,jj) = (y_k(ii,1)*x_v(jj,1)-y_v(jj,1))/(2*pi*r_v(jj,1)^2);
end
end
end
function xgamma_mat = solve(obj)
xy_bb = get_bb(obj);
b_mat = xy_bb(:,2);
A_mat = get_AA(obj);
gamma = A_mat\b_mat;
objMesh = mesh(obj.nCells,obj.f_max);
deltax = get_dx(objMesh);
gamma = gamma ./ deltax;
%for ii = 1 : size(gamma,1)
% deltac = get_dc(objMesh,ii);
% gamma(ii,1) = gamma(ii,1) ./ deltac;
%end
xy_m = mesh_vortex_points(objMesh);
x_m = xy_m(:,1);
xgamma_mat = [x_m,gamma];
end
function solve_plot(obj)
xgamma_mat = solve(obj);
figure(1);
alpha_deg = obj.alpha * 180 / pi;
objExt = vortex_ext(obj.f_max,alpha_deg);
A_ext = ext_matrix(objExt);
x_ext = A_ext(:,1);
y_ext = A_ext(:,2);
plot(x_ext,y_ext); hold on;
x_num = xgamma_mat(:,1);
y_num = xgamma_mat(:,2);
scatter(x_num,y_num); hold on;
end
end
end
VLMalphaList.m
% VLMalphaList.m
classdef VLMalphaList
properties
f_max
alpha_deg_list
nCellsList
end
methods
function obj = VLMalphaList(f_max,alpha_deg_list,nCellsList)
obj.f_max = f_max;
obj.alpha_deg_list = alpha_deg_list;
obj.nCellsList = nCellsList;
end
function solve_plot(obj)
for ii = 1 : size(obj.alpha_deg_list,1)
objVLMList = VLMList(obj.f_max,obj.alpha_deg_list(ii,1),obj.nCellsList);
solve_plot(objVLMList);
end
end
end
end
VLMList.m
% VLMList.m
classdef VLMList
properties
f_max
alpha
nCellsList
end
methods
function obj = VLMList(f_max,alpha_deg,nCellsList)
obj.f_max = f_max;
obj.alpha = alpha_deg / 180 * pi;
obj.nCellsList = nCellsList;
end
function solve_plot(obj)
alpha_deg = obj.alpha * 180 / pi;
objExt = vortex_ext(obj.f_max,alpha_deg);
xy_ext = ext_matrix(objExt);
x_ext = xy_ext(:,1);
y_ext = xy_ext(:,2);
figure(1);
plt0 = plot(x_ext,y_ext); hold on;
plt0.Color = 'black';
plt0.LineWidth = 1.5;
nCellsListLength = size(obj.nCellsList,1);
for ii = 1 : size(obj.nCellsList,1)
objVLM = VLM(obj.f_max,alpha_deg,obj.nCellsList(ii,1));
xgamma_mat = solve(objVLM);
x_num = xgamma_mat(:,1);
y_num = xgamma_mat(:,2);
plt1=plot(x_num,y_num); hold on;
plt1.Color = [0,0,0,(ii/nCellsListLength)];
plt1.LineStyle = '--';
plt1.LineWidth = 1.0;
sc1=scatter(x_num,y_num); hold on;
sc1.Marker = 'x';
sc1.MarkerEdgeColor = [0,0,0];
sc1.MarkerEdgeAlpha = ii/nCellsListLength;
sc1.SizeData = 400/obj.nCellsList(ii,1);
end
plt1=plot(x_num,y_num); hold on;
plt1.Color = [0,0,0,(ii/nCellsListLength)];
plt1.LineStyle = '--';
plt1.LineWidth = 1.5;
plt1.Marker = 'x';
xL1 = xlabel('x');
xL1.FontSize = 17;
yL1 = ylabel('\gamma(x)');
yL1.FontSize = 17;
leg1 = legend([plt0,plt1]);
leg1.String = ["Ext","Num"];
leg1.FontSize = 17;
leg1.Location = 'best';
leg1.Color = 'none';
titleStr = "$f_{\max} = $ " + num2str(obj.f_max) + "\,, ";
titleStr = titleStr + "$\alpha = $ " + num2str(alpha_deg) + "$\,^\circ$";
tit1 = title(titleStr);
tit1.FontSize = 17;
tit1.Interpreter = 'latex';
saveStr = "./figs/fg_";
saveStr = saveStr + "f" + num2str(obj.f_max);
saveStr = saveStr + "a" + num2str(alpha_deg);
savePng = saveStr + ".png";
saveas(gcf,savePng);
saveEps = saveStr + ".eps";
saveas(gcf,saveEps,'epsc2');
close all;
end
end
end
vortex_ext.m
% vortex_ext.m
classdef vortex_ext
properties
f_max=0.08;
alpha_deg = 5;
end
methods
function obj = vortex_ext(f_max,alpha_deg)
obj.f_max = f_max;
obj.alpha_deg = alpha_deg;
end
function A_ext = ext_matrix(obj)
c=1;U=1;
alpha = obj.alpha_deg / 180 * pi;
theta=0:0.0001:pi;
x=zeros(1,size(theta,2));
gamma=zeros(1,size(x,2));
x_min=0.01;
for i = 1:size(x,2)
x(1,i)=get_x(theta(1,i));
gamma(1,i)=get_gamma(alpha,theta(1,i),obj.f_max);
end
for i = 1:size(x,2)
if(x(1,i)>=x_min)
index1 = i;
break;
end
end
for i = 1:size(x,2)
x_max = 1-x_min;
if(x(1,i)>=x_max)
index2 = i;
break;
end
end
x=x(1,(index1:index2));
gamma=gamma(1,(index1:index2));
x_ext = x';
gamma_ext = gamma';
A_ext = [x_ext,gamma_ext];
function x=get_x(theta)
c = 1;
x = (c/2)*(1-cos(theta));
end
function gamma=get_gamma(alpha,theta,f_max)
c = 1;
U = 1;
gamma1 = 2*U*alpha*((1+cos(theta))/sin(theta));
gamma2 = 8*U*(f_max/c)*sin(theta);
gamma = gamma1 + gamma2;
end
end
end
end
main.m
clearvars; clc; % main.m
alphaDegList = [0;4;8];
nCellsList = [200;500;800];
objVLMSalphaList = VLMSalphaList(alphaDegList,nCellsList);
plot_solve(objVLMSalphaList);
meshS.m
% meshS.m
classdef meshS
properties
mesh_xyy
nCells
end
methods
function obj = meshS(nCells)
obj.nCells = nCells;
xyy_org = [0,0,0;0.0125,0.0244,-0.0143;0.025,0.0339,-0.0195;0.05,0.0473,-0.0249;...
0.075,0.0576,-0.0274;0.1,0.0659,-0.0286;0.15,0.0789,-0.0288;0.2,0.0880,-0.0274;...
0.25,0.0941,-0.0250;0.3,0.0976,-0.0226;0.4,0.0980,-0.0180;0.5,0.0919,-0.0140;...
0.6,0.0814,-0.0100;0.7,0.0669,-0.0065;0.8,0.0489,-0.0039;0.9,0.0271,-0.0022;...
0.95,0.0147,-0.0016;1,0.0013,-0.0013];
x_org = xyy_org(:,1);
y1_org = xyy_org(:,2);
y2_org = xyy_org(:,3);
dx = 1/obj.nCells;
x_int = 0.01:dx:1;
x_int = [0,x_int];
x_int = x_int';
int_method = 'spline'; %linear'spline'
y1_int=interp1(x_org,y1_org,x_int,int_method);
y2_int=interp1(x_org,y2_org,x_int,int_method);
obj.mesh_xyy = [x_int,y1_int,y2_int];
end
function upper = get_upper(obj)
x_upper = obj.mesh_xyy(:,1);
y_upper = obj.mesh_xyy(:,2);
upper = [x_upper,y_upper];
end
function lower = get_lower(obj)
x_lower = obj.mesh_xyy(:,1);
y_lower = obj.mesh_xyy(:,3);
lower = [x_lower,y_lower];
end
function middle = get_middle(obj)
upper = get_upper(obj);
lower = get_lower(obj);
y_upper = upper(:,2);
y_lower = lower(:,2);
x_middle = obj.mesh_xyy(:,1);
y_middle = zeros(size(y_upper,1),1);
for ii = 1 : size(y_middle,1)
y_middle(ii,1) = (y_upper(ii,1) + y_lower(ii,1))/2;
end
middle = [x_middle,y_middle];
end
function plot_foil(obj)
upper = get_upper(obj);
lower = get_lower(obj);
middle = get_middle(obj);
xllim = -0.1;
xrlim = 1.1;
yllim = -0.3;
yrlim = 0.4;
figure(1);
set(gca,'xlim',[xllim,xrlim]);
hold on;
set(gca,'ylim',[yllim,yrlim]);
hold on;
x_upper = upper(:,1);
y_upper = upper(:,2);
plt_upper = plot(x_upper,y_upper);
plt_upper.LineWidth = 1.5;
plt_upper.Color = 'blue';
plt_upper.Marker = 'x';
hold on;
x_lower = lower(:,1);
y_lower = lower(:,2);
plt_lower = plot(x_lower,y_lower);
plt_lower.LineWidth = 1.5;
plt_lower.Color = 'blue';
plt_lower.Marker = 'x';
hold on;
x_middle = middle(:,1);
y_middle = middle(:,2);
plt_middle = plot(x_middle,y_middle);
plt_middle.LineWidth = 1.5;
plt_middle.Color = 'red';
plt_middle.Marker = 'x';
hold on;
end
function thick = get_thick(obj)
x_thick = obj.mesh_xyy(:,1);
y_thick = zeros(size(x_thick,1),1);
for ii = 1 : size(y_thick,1)
y_thick(ii,1) = obj.mesh_xyy(ii,2) - obj.mesh_xyy(ii,3);
end
thick = [x_thick,y_thick];
end
function xy_p = mesh_points(obj)
middle = get_middle(obj);
xy_p = middle;
end
function slope = mesh_slope(obj)
xy_p = mesh_points(obj);
x_p = xy_p(:,1);
y_p = xy_p(:,2);
y_k = zeros(size(y_p,1)-1,1);
x_k = zeros(size(y_k,1),1);
for ii = 1 : size(x_k,1)
y_k(ii,1) = (y_p(ii+1,1)-y_p(ii,1))/(x_p(ii+1,1)-x_p(ii,1));
x_k(ii,1) = (x_p(ii+1,1)+x_p(ii,1))/2;
end
slope = [x_k,y_k];
end
function linear = mesh_linear_points(obj,fraction)
xy_p = mesh_points(obj);
xy_k = mesh_slope(obj);
y_k = xy_k(:,2);
x_p = xy_p(:,1);
y_p = xy_p(:,2);
x_line = zeros(size(y_k,1),1);
y_line = zeros(size(x_line,1),1);
for ii = 1 : size(y_line,1)
delta_x = x_p(ii+1,1) - x_p(ii,1);
delta_x = delta_x*fraction;
x_line(ii,1) = x_p(ii,1) + delta_x;
y_line(ii,1) = y_p(ii,1) + delta_x*y_k(ii,1);
end
linear = [x_line,y_line];
end
function xy_mat = mesh_vortex_points(obj)
fraction = 0.25;
xy_mat = mesh_linear_points(obj,fraction);
end
function xy_mat = mesh_middle_points(obj)
fraction = 0.5;
xy_mat = mesh_linear_points(obj,fraction);
end
function xy_mat = mesh_control_points(obj)
fraction = 0.75;
xy_mat = mesh_linear_points(obj,fraction);
end
function xyr_mat = mesh_vortex_xyr(obj,xp,yp)
xy_v = mesh_vortex_points(obj);
x_v = xy_v(:,1);
y_v = xy_v(:,2);
x_mat = zeros(size(y_v,1),1);
y_mat = zeros(size(x_mat,1),1);
r_mat = zeros(size(y_mat,1),1);
for ii = 1 : size(r_mat,1)
xx = xp - x_v(ii,1);
yy = yp - y_v(ii,1);
r_mat(ii,1) = sqrt(xx^2+yy^2);
x_mat(ii,1) = yy;
y_mat(ii,1) = - xx;
end
xyr_mat = [x_mat,y_mat,r_mat];
end
function deltac = get_deltac(obj)
xy_p = mesh_points(obj);
x_p = xy_p(:,1);
y_p = xy_p(:,2);
deltac = zeros(size(y_p,1)-1,1);
for ii = 1 : size(deltac,1)
dx = x_p(ii+1,1) - x_p(ii,1);
dy = y_p(ii+1,1) - y_p(ii,1);
deltac(ii,1) = sqrt(dx^2+dget_qy^2);
end
end
function qq = get_q(obj)
U = 1;
thick = get_thick(obj);
y_t = thick(:,2);
xyr_tan = get_xyr_tan(obj);
x_tan = xyr_tan(:,1);
r_tan = xyr_tan(:,3);
qq = zeros(size(y_t,1)-1,1);
for ii = 1 : size(qq,1)
dt = y_t(ii+1,1) - y_t(ii,1);
dx = x_tan(ii,1);
delta_x = r_tan(ii,1);
qq(ii,1) = U * (dt/dx) * delta_x;
end
end
function xyr_mat = mesh_q_xyr(obj,xp,yp)
xy_v = mesh_vortex_points(obj);
x_v = xy_v(:,1);
y_v = xy_v(:,2);
x_mat = zeros(size(y_v,1),1);
y_mat = zeros(size(x_mat,1),1);
r_mat = zeros(size(y_mat,1),1);
for ii = 1 : size(r_mat,1)
x_mat(ii,1) = xp - x_v(ii,1);
y_mat(ii,1) = yp - y_v(ii,1);
r_mat(ii,1) = sqrt(x_mat(ii,1)^2+y_mat(ii,1)^2);
end
xyr_mat = [x_mat,y_mat,r_mat];
end
function xyr_tan = get_xyr_tan(obj)
xy_p = mesh_points(obj);
x_p = xy_p(:,1);
y_p = xy_p(:,2);
x_tan = zeros(size(y_p,1)-1,1);
y_tan = zeros(size(x_tan,1),1);
r_tan = zeros(size(y_tan,1),1);
for ii = 1 : size(r_tan,1)
xx = x_p(ii+1,1) - x_p(ii,1);
yy = y_p(ii+1,1) - y_p(ii,1);
r_tan(ii,1) = sqrt(xx^2+yy^2);
x_tan(ii,1) = xx;
y_tan(ii,1) = yy;
end
xyr_tan = [x_tan,y_tan,r_tan];
end
function xyr_tan = get_xyr_tan_upper(obj)
xy_p = get_upper(obj);
x_p = xy_p(:,1);
y_p = xy_p(:,2);
x_tan = zeros(size(y_p,1)-1,1);
y_tan = zeros(size(x_tan,1),1);
r_tan = zeros(size(y_tan,1),1);
for ii = 1 : size(r_tan,1)
xx = x_p(ii+1,1) - x_p(ii,1);
yy = y_p(ii+1,1) - y_p(ii,1);
r_tan(ii,1) = sqrt(xx^2+yy^2);
x_tan(ii,1) = xx;
y_tan(ii,1) = yy;
end
xyr_tan = [x_tan,y_tan,r_tan];
end
function xyr_tan = get_xyr_tan_lower(obj)
xy_p = get_lower(obj);
x_p = xy_p(:,1);
y_p = xy_p(:,2);
x_tan = zeros(size(y_p,1)-1,1);
y_tan = zeros(size(x_tan,1),1);
r_tan = zeros(size(y_tan,1),1);
for ii = 1 : size(r_tan,1)
xx = x_p(ii+1,1) - x_p(ii,1);
yy = y_p(ii+1,1) - y_p(ii,1);
r_tan(ii,1) = sqrt(xx^2+yy^2);
x_tan(ii,1) = xx;
y_tan(ii,1) = yy;
end
xyr_tan = [x_tan,y_tan,r_tan];
end
end
end
nacaExp.m
% nacaExp.m
classdef nacaExp
properties
xy_u
xy_l
end
methods
function obj = nacaExp()
xy_upper = [0.0020638925306921407, -2.145504423413796;...
0.00930955082794406, -2.1885817532803777;...
0.016400968046389852, -2.145592561173114;...
0.03061685414302555, -2.0780570030958394;...
0.04842802633850038, -2.0166910881708113;...
0.07518885351137161, -1.9492326506329396;...
0.10019059790452475, -1.9002060220124055;...
0.1252474283972516, -1.8819174369539593;...
0.14851212445051612, -1.863617834675598;...
0.17532803772296093, -1.8268974406998142;...
0.20038486821568782, -1.8086088556413678;...
0.2236605814888671, -1.796456862075424;...
0.2469032431023022, -1.765862042372228;...
0.2737411908145765, -1.7414368658212789;...
0.30236025574639824, -1.710875097777827;...
0.3470534445338064, -1.64967444115152;...
0.3953088677603094, -1.5762006015402075;...
0.4490067976246873, -1.5396454658631444;...
0.5008134380703706, -1.447750834554409;...
0.5490908957367031, -1.3865722123679314;...
0.6010077083815337, -1.3561536681833708;...
0.6510993349271578, -1.307281280641643;...
0.7029500442525001, -1.239977084182578;...
0.7548227880176716, -1.184968105148347;...
0.8031222801238336, -1.1360847003867045;...
0.8478264861311563, -1.0810316524728152;...
0.8996771954564985, -1.0137274560137497;...
0.9443593670239918, -0.9463791906750254;...
0.9818950352734657, -0.8912820738814768;...
0.9980352624485405, -0.8975288375731272;...
];
xy_lower = [-0.002008806431118504, -0.8729384027234572;...
0.006070488173014407, -0.38118479182963005;...
0.01311783651180129, -0.31360516487269674;...
0.01680126036995823, -0.3689556777242835;...
0.02956654584448827, -0.49198497251203666;...
0.05129250351632936, -0.6150693533993636;...
0.07827367508749511, -0.6705630901098423;...
0.09990047778010361, -0.7383189925854117;...
0.15370857984362893, -0.7632399440325233;...
0.20034814414930538, -0.7881168265999765;...
0.25052790845424733, -0.788425308757589;...
0.3007297071990187, -0.8010290083400355;...
0.34910631984458373, -0.7951788645653162;...
0.3975270013698077, -0.8139191556402663;...
0.45309418621305086, -0.8204082981700402;...
0.5014707988586159, -0.8145581543953204;...
0.5462741598451714, -0.8148335848931887;...
0.6520321262132713, -0.8277788182929924;...
0.7542058237024468, -0.8407020172529664;...
0.8527842350927833, -0.8474555730606941;...
0.9029860338375546, -0.860059272643142;...
0.9442161431651004, -0.8664602774135979;...
0.9728902941964959, -0.8666365529322335;...
0.9962541452289928, -0.9036654290656299;...
];
y_upper = xy_upper(:,2);
y_upper = y_upper + ones(size(y_upper,1),1) .* 0.75;
y_upper = y_upper ./3;
y_upper = y_upper .* 4.5;
x_upper = xy_upper(:,1);
xy_upper = [x_upper,y_upper];
y_lower = xy_lower(:,2);
y_lower = y_lower + ones(size(y_lower,1),1) .* 0.75;
y_lower = y_lower ./3;
y_lower = y_lower .* 4.5;
x_lower = xy_lower(:,1);
xy_lower = [x_lower,y_lower];
obj.xy_u = xy_upper;
obj.xy_l = xy_lower;
end
end
end
VLMS.m
% VLMS.m
classdef VLMS
properties
alpha
ifUpper
nCells
end
methods
function obj = VLMS(alpha_deg,nCells)
obj.alpha = alpha_deg / 180 * pi;
obj.nCells = nCells;
end
function y_bb = get_bb(obj)
U = 1;
objMesh = meshS(obj.nCells);
xy_k = mesh_slope(objMesh);
y_k = xy_k(:,2);
xy_c = mesh_control_points(objMesh);
x_c = xy_c(:,1);
y_c = xy_c(:,2);
y_bb = zeros(size(y_k,1),1);
Ux = U * cos(obj.alpha);
Uy = U * sin(obj.alpha);
for ii = 1 : size(y_bb,1)
Qx = 0.0;
Qy = 0.0;
xyr_q = mesh_q_xyr(objMesh,x_c(ii,1),y_c(ii,1));
qq = get_q(objMesh);
x_q = xyr_q(:,1);
y_q = xyr_q(:,2);
r_q = xyr_q(:,3);
for jj = 1 : size(r_q,1)
Qx = Qx + (qq(jj,1)*x_q(jj,1))/(2*pi*r_q(jj,1)^2);
Qy = Qy + (qq(jj,1)*y_q(jj,1))/(2*pi*r_q(jj,1)^2);
end
y_bb(ii,1) = (Uy + Qy) - y_k(ii,1)*(Ux+Qx);
end
end
function A_mat = get_AA(obj)
objMesh = meshS(obj.nCells);
xy_c = mesh_control_points(objMesh);
xy_k = mesh_slope(objMesh);
y_k = xy_k(:,2);
x_c = xy_c(:,1);
y_c = xy_c(:,2);
A_mat = zeros(size(y_c,1),size(y_c,1));
for ii = 1 : size(A_mat,1)
xp = x_c(ii,1);
yp = y_c(ii,1);
xyr_v = mesh_vortex_xyr(objMesh,xp,yp);
x_v = xyr_v(:,1);
y_v = xyr_v(:,2);
r_v = xyr_v(:,3);
for jj = 1 : size(A_mat,2)
A_mat(ii,jj) = (y_k(ii,1)*x_v(jj,1)-y_v(jj,1))/(2*pi*r_v(jj,1)^2);
end
end
end
function gamma = get_gamma(obj)
b_mat = get_bb(obj);
A_mat = get_AA(obj);
gamma = A_mat\b_mat;
end
function vt0 = get_vt0(obj)
U = 1;
objMesh = meshS(obj.nCells);
xyr_vt0 = get_xyr_tan(objMesh);
x_vt0 = xyr_vt0(:,1);
y_vt0 = xyr_vt0(:,2);
r_vt0 = xyr_vt0(:,3);
vt0 = zeros(size(r_vt0,1),1);
Ux = U * cos(obj.alpha);
Uy = U * sin(obj.alpha);
for ii = 1 : size(vt0,1)
vt0(ii,1) = (Ux*x_vt0(ii,1)+Uy*y_vt0(ii,1))/r_vt0(ii,1);
end
end
function xy_vt1_vortex = get_vt1_vortex(obj)
gamma = get_gamma(obj);
objMesh = meshS(obj.nCells);
xy_c = mesh_control_points(objMesh);
x_c = xy_c(:,1);
y_c = xy_c(:,2);
x_vt1_vortex = zeros(size(y_c,1),1);
y_vt1_vortex = zeros(size(x_vt1_vortex,1),1);
for ii = 1 : size(y_vt1_vortex,1)
xp = x_c(ii,1);
yp = y_c(ii,1);
xyr_v = mesh_vortex_xyr(objMesh,xp,yp);
x_v = xyr_v(:,1);
y_v = xyr_v(:,2);
r_v = xyr_v(:,3);
x_vt1 = 0.0;
y_vt1 = 0.0;
for jj = 1 : size(r_v,1)
if (jj ~= ii)
x_vt1_jj = (gamma(jj,1)*x_v(jj,1))/(2*pi*r_v(jj,1)^2);
x_vt1 = x_vt1 + x_vt1_jj;
y_vt1_jj = (gamma(jj,1)*y_v(jj,1))/(2*pi*r_v(jj,1)^2);
y_vt1 = y_vt1 + y_vt1_jj;
end
end
x_vt1_vortex(ii,1) = x_vt1;
y_vt1_vortex(ii,1) = y_vt1;
end
xy_vt1_vortex = [x_vt1_vortex,y_vt1_vortex];
end
function xy_vt1_source = get_vt1_source(obj)
objMesh = meshS(obj.nCells);
xy_c = mesh_control_points(objMesh);
x_c = xy_c(:,1);
y_c = xy_c(:,2);
x_vt1_source = zeros(size(y_c,1),1);
y_vt1_source = zeros(size(x_vt1_source,1),1);
for ii = 1 : size(y_vt1_source,1)
xp = x_c(ii,1);
yp = y_c(ii,1);
xyr_q = mesh_q_xyr(objMesh,xp,yp);
qq = get_q(objMesh);
x_q = xyr_q(:,1);
y_q = xyr_q(:,2);
r_q = xyr_q(:,3);
Qx = 0.0;
Qy = 0.0;
for jj = 1 : size(r_q,1)
if (jj ~= ii)
Qx_jj = (qq(jj,1)*x_q(jj,1))/(2*pi*r_q(jj,1)^2);
Qx = Qx + Qx_jj;
Qy_jj = (qq(jj,1)*y_q(jj,1))/(2*pi*r_q(jj,1)^2);
Qy = Qy + Qy_jj;
end
end
x_vt1_source(ii,1) = Qx;
y_vt1_source(ii,1) = Qy;
end
xy_vt1_source = [x_vt1_source,y_vt1_source];
end
function vt1 = get_vt1(obj)
objMesh = meshS(obj.nCells);
xyr_tan = get_xyr_tan(objMesh);
x_tan = xyr_tan(:,1);
y_tan = xyr_tan(:,2);
r_tan = xyr_tan(:,3);
xy_vt1_source = get_vt1_source(obj);
x_source = xy_vt1_source(:,1);
y_source = xy_vt1_source(:,2);
xy_vt1_vortex = get_vt1_vortex(obj);
x_vortex = xy_vt1_vortex(:,1);
y_vortex = xy_vt1_vortex(:,2);
x_vt1 = x_source + x_vortex;
y_vt1 = y_source + y_vortex;
vt1 = zeros(size(y_vt1,1),1);
for ii = 1 : size(vt1,1)
vt1(ii,1) = (x_vt1(ii,1)*x_tan(ii,1)+y_vt1(ii,1)*y_tan(ii,1))/r_tan(ii,1);
end
end
function vt2 = get_vt2(obj)
gamma = get_gamma(obj);
vt2 = zeros(size(gamma,1),1);
objMesh = meshS(obj.nCells);
xy_p = mesh_points(objMesh);
x_p = xy_p(:,1);
for ii = 1 : size(vt2,1)
dx = x_p(ii+1,1) - x_p(ii,1);
vt2(ii,1) = 0.5 * gamma(ii,1) / dx;
end
end
function Cp = get_Cp(obj,fraction)
U = 1;
vt0 = get_vt0(obj);
vt1 = get_vt1(obj);
vt2 = get_vt2(obj);
Cp = zeros(size(vt2,1),1);
for ii = 1 : size(Cp,1)
sumV = vt0(ii,1) + vt1(ii,1) + fraction*vt2(ii,1);
Cp(ii,1) = 1 - (sumV/U)^2;
end
end
function Cpu = get_Cpu(obj)
Cpu = get_Cp(obj,1.0);
end
function Cpl = get_Cpl(obj)
Cpl = get_Cp(obj,-1.0);
end
function plot_cp(obj)
Cpu = get_Cpu(obj);
Cpl = get_Cpl(obj);
objMesh = meshS(obj.nCells);
xy_c = mesh_control_points(objMesh);
x_c = xy_c(:,1);
figure(1);
plt_u = plot(x_c,Cpu);
plt_u.LineWidth = 1.5;
plt_u.Color = 'blue';
hold on;
plt_l = plot(x_c,Cpl);
plt_l.LineWidth = 1.5;
plt_l.Color = 'red';
hold on;
end
function Gamma = get_Gamma(obj)
gamma = get_gamma(obj);
objMesh = meshS(obj.nCells);
xyr_tan = get_xyr_tan(objMesh);
r_tan = xyr_tan(:,3);
Gamma = zeros(size(gamma,1),1);
for ii = 1 : size(Gamma,1)
deltax = r_tan(ii,1);
Gamma(ii,1) = gamma(ii,1) / deltax;
end
end
function xgamma = get_xgamma(obj)
gamma = get_Gamma(obj);
objMesh = meshS(obj.nCells);
xy_mat = mesh_vortex_points(objMesh);
x_mat = xy_mat(:,1);
xgamma = [x_mat,gamma];
end
function Source = get_Source(obj)
objMesh = meshS(obj.nCells);
source = get_q(objMesh);
xyr_tan = get_xyr_tan(objMesh);
r_tan = xyr_tan(:,3);
Source = zeros(size(source,1),1);
for ii = 1 : size(Source,1)
deltax = r_tan(ii,1);
Source(ii,1) = source(ii,1) / deltax;
end
end
function plot_cp_ifUpper(obj)
objExp = nacaExp();
x_extu = objExp.xy_u(:,1);
y_extu = objExp.xy_u(:,2);
x_extl = objExp.xy_l(:,1);
y_extl = objExp.xy_l(:,2);
Cpu = get_Cpu_ifUpper(obj);
Cpl = get_Cpl_ifUpper(obj);
%Cpu = get_Cpu(obj);
%Cpl = get_Cpl(obj);
objMesh = meshS(obj.nCells);
xy_c = mesh_control_points(objMesh);
x_c = xy_c(:,1);
figure(1);
set(gca,'xlim',[0,1]);
hold on;
set(gca,'ylim',[-3.5,1.5]);
hold on;
plt_u = plot(x_c,Cpu);
plt_u.LineWidth = 1.5;
plt_u.Color = 'blue';
hold on;
plt_l = plot(x_c,Cpl);
plt_l.LineWidth = 1.5;
plt_l.Color = 'red';
hold on;
sc1 = scatter(x_extu,y_extu);
sc1.Marker = 'x';
sc1.MarkerEdgeColor = 'blue';
hold on;
sc2 = scatter(x_extl,y_extl);
sc2.Marker = 'x';
sc2.MarkerEdgeColor = 'red';
hold on;
end
function Cpu = get_Cpu_ifUpper(obj)
obj.ifUpper = 1;
Cpu = get_Cp_ifUpper(obj);
end
function Cpl = get_Cpl_ifUpper(obj)
obj.ifUpper = 0;
Cpl = get_Cp_ifUpper(obj);
end
function Cp = get_Cp_ifUpper(obj)
if(obj.ifUpper==1)
fraction = 1.0;
elseif(obj.ifUpper==0)
fraction = -1.0;
end
U = 1;
vt0 = get_vt0_ifUpper(obj);
vt1 = get_vt1_ifUpper(obj);
vt2 = get_vt2(obj);
Cp = zeros(size(vt2,1),1);
for ii = 1 : size(Cp,1)
sumV = vt0(ii,1) + vt1(ii,1) + fraction*vt2(ii,1);
Cp(ii,1) = 1 - (sumV/U)^2;
end
end
function vt0 = get_vt0_ifUpper(obj)
U = 1;
objMesh = meshS(obj.nCells);
if(obj.ifUpper==1)
xyr_vt0 = get_xyr_tan_upper(objMesh);
elseif(obj.ifUpper==0)
xyr_vt0 = get_xyr_tan_lower(objMesh);
end
x_vt0 = xyr_vt0(:,1);
y_vt0 = xyr_vt0(:,2);
r_vt0 = xyr_vt0(:,3);
vt0 = zeros(size(r_vt0,1),1);
Ux = U * cos(obj.alpha);
Uy = U * sin(obj.alpha);
for ii = 1 : size(vt0,1)
vt0(ii,1) = (Ux*x_vt0(ii,1)+Uy*y_vt0(ii,1))/r_vt0(ii,1);
end
end
function vt1 = get_vt1_ifUpper(obj)
objMesh = meshS(obj.nCells);
if(obj.ifUpper==1)
xyr_tan = get_xyr_tan_upper(objMesh);
elseif(obj.ifUpper==0)
xyr_tan = get_xyr_tan_lower(objMesh);
end
x_tan = xyr_tan(:,1);
y_tan = xyr_tan(:,2);
r_tan = xyr_tan(:,3);
xy_vt1_source = get_vt1_source(obj);
x_source = xy_vt1_source(:,1);
y_source = xy_vt1_source(:,2);
xy_vt1_vortex = get_vt1_vortex(obj);
x_vortex = xy_vt1_vortex(:,1);
y_vortex = xy_vt1_vortex(:,2);
x_vt1 = x_source + x_vortex;
y_vt1 = y_source + y_vortex;
vt1 = zeros(size(y_vt1,1),1);
for ii = 1 : size(vt1,1)
vt1(ii,1) = (x_vt1(ii,1)*x_tan(ii,1)+y_vt1(ii,1)*y_tan(ii,1))/r_tan(ii,1);
end
end
end
end
VLMSalphaList.m
% VLMSalphaList.m
classdef VLMSalphaList
properties
alphaDegList
nCellsList
end
methods
function obj = VLMSalphaList(alphaDegList,nCellsList)
obj.alphaDegList = alphaDegList;
obj.nCellsList = nCellsList;
end
function plot_solve(obj)
for ii = 1 : size(obj.alphaDegList,1)
alpha_deg = obj.alphaDegList(ii,1);
objVLMSnCellsList = VLMSnCellsList(alpha_deg,obj.nCellsList);
plot_gamma(objVLMSnCellsList);
plot_cpul(objVLMSnCellsList);
end
end
end
end
VLMSnCellsList.m
% VLMSnCellsList.m
classdef VLMSnCellsList
properties
alpha_deg
nCellsList
end
methods
function obj = VLMSnCellsList(alpha_deg,nCellsList)
obj.alpha_deg = alpha_deg;
obj.nCellsList = nCellsList;
end
function plot_gamma(obj)
figure(1);
lineStyle = [":";"--";"-"];
plts = zeros(size(obj.nCellsList,1),1);
for ii = 1 : size(obj.nCellsList,1)
objVLMS = VLMS(obj.alpha_deg,obj.nCellsList(ii,1));
xgamma = get_xgamma(objVLMS);
x = xgamma(:,1);
gamma = xgamma(:,2);
plt1=plot(x,gamma);
plt1.Color = 'black';
plt1.LineWidth = 1.5;
plt1.LineStyle = lineStyle(ii,1);
hold on;
plts(ii,1) = plt1;
end
leg1 = legend(plts);
leg1.FontSize = 17;
leg1.Location = 'best';
leg1.Color = 'none';
leg1.String = ["N=200","N=500","N=800"];
xL1 = xlabel('x');
xL1.FontSize = 17;
yL1 = ylabel('\gamma(x)');
yL1.FontSize = 17;
titleStr = "$\alpha = $ " + num2str(obj.alpha_deg) + "$\,^\circ$";
tit1 = title(titleStr);
tit1.FontSize = 17;
tit1.Interpreter = 'latex';
saveStr = "./figs/gx_";
saveStr = saveStr + "a" + num2str(obj.alpha_deg);
savePng = saveStr + ".png";
saveas(gcf,savePng);
saveEps = saveStr + ".eps";
saveas(gcf,saveEps,'epsc2');
close all;
end
function plot_cpul(obj)
figure(1);
objExp = nacaExp();
x_extu = objExp.xy_u(:,1);
y_extu = objExp.xy_u(:,2);
x_extl = objExp.xy_l(:,1);
y_extl = objExp.xy_l(:,2);
set(gca,'xlim',[0,1]);
hold on;
set(gca,'ylim',[-3.5,1.5]);
hold on;
if (obj.alpha_deg == 8)
sc1 = scatter(x_extu,y_extu);
sc1.Marker = 'x';
sc1.MarkerEdgeColor = 'blue';
sc1.SizeData = 80;
hold on;
sc2 = scatter(x_extl,y_extl);
sc2.Marker = 'x';
sc2.MarkerEdgeColor = 'red';
sc2.SizeData = 80;
hold on;
end
lineStyle = [":";"--";"-"];
plts = zeros(size(obj.nCellsList,1),2);
for ii = 1 : size(obj.nCellsList,1)
objVLMS = VLMS(obj.alpha_deg,obj.nCellsList(ii,1));
Cpu = get_Cpu_ifUpper(objVLMS);
Cpl = get_Cpl_ifUpper(objVLMS);
objMesh = meshS(obj.nCellsList(ii,1));
xy_c = mesh_control_points(objMesh);
x_c = xy_c(:,1);
plt1 = plot(x_c,Cpu);
plt1.LineWidth = 1.5;
plt1.Color = 'blue';
plt1.LineStyle = lineStyle(ii,1);
hold on;
plt2 = plot(x_c,Cpl);
plt2.LineWidth = 1.5;
plt2.Color = 'red';
plt2.LineStyle = lineStyle(ii,1);
hold on;
plts(ii,1) = plt1;
plts(ii,2) = plt2;
end
if (obj.alpha_deg == 8)
leg1 = legend([sc2,plts(3,1)]);
leg1.String = ["Exp.lower","Num.upper"];
else
leg1 = legend([plts(1,1),plts(2,2),plts(3,1)]);
leg1.String = ["N=200.upper","N=500.lower","N=800.upper"];
end
leg1.FontSize = 17;
leg1.Location = 'best';
leg1.Color = 'none';
xL1 = xlabel('x');
xL1.FontSize = 17;
yL1 = ylabel('Cp');
yL1.FontSize = 17;
titleStr = "$\alpha = $ " + num2str(obj.alpha_deg) + "$\,^\circ$";
tit1 = title(titleStr);
tit1.FontSize = 17;
tit1.Interpreter = 'latex';
saveStr = "./figs/cp_";
saveStr = saveStr + "a" + num2str(obj.alpha_deg);
savePng = saveStr + ".png";
saveas(gcf,savePng);
saveEps = saveStr + ".eps";
saveas(gcf,saveEps,'epsc2');
close all;
end
end
end
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
