G11-NAOE9101MSP-船舶推进大作业

最新推荐文章于 2025-07-05 11:15:07 发布

jmsyh

最新推荐文章于 2025-07-05 11:15:07 发布

阅读量160

点赞数

CC 4.0 BY-SA版权

文章标签： javascript python 开发语言

本文链接：https://blog.youkuaiyun.com/joshua_shi_t/article/details/128989093

该系列MATLAB代码实现了一种变稳方法（VortexLatticeMethod,VLM），用于计算机翼的气动特性。主要函数包括`VLM.m`、`mesh.m`、`VLMalphaList.m`等，涉及网格生成、涡旋分布计算、气动力计算等步骤。代码首先定义了不同角度攻角（alpha_deg）和网格数（nCells）的组合，然后对每个组合进行VLM计算并绘制结果，展示了不同攻角和网格数量对机翼升力和阻力的影响。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

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