给出一份基于matlab的计算绕圆柱无环量流动绘制流线的程序代码

以下是一个简单的基于MATLAB的计算绕圆柱无环量流动绘制流线的程序代码：

clear all;
clc;

% 定义常数
R = 1; % 圆柱半径
U = 1; % 入口速度
N = 50; % 流线数量
M = 100; % 离散点数量

% 定义圆柱上的离散点坐标
theta = linspace(0, 2*pi, M);
x = R*cos(theta);
y = R*sin(theta);

% 定义网格坐标
x_min = -2;
x_max = 2;
y_min = -2;
y_max = 2;
dx = 0.1;
dy = 0.1;
[x_grid, y_grid] = meshgrid(x_min:dx:x_max, y_min:dy:y_max);

% 计算速度势函数和流函数
phi = U*x_grid;
psi = U*y_grid;

for i = 1:N % 针对每一条流线
    % 定义起点坐标
    x_start = -2;
    y_start = y_max/(N+1)*i;
    
    % 定义流线终点坐标
    x_end = x_max;
    y_end = y_start;
    
    % 计算流线轨迹
    x_streamline = linspace(x_start, x_end, M);
    y_streamline = linspace(y_start, y_end, M);
    for j = 2:M
        x_current = x_streamline(j);
        y_current = y_streamline(j);
        
        % 计算速度势函数和流函数
        phi_current = U*x_current;
        psi_current = U*y_current;
        for k = 1:M
            x_p = x(k);
            y_p = y(k);
            r = sqrt((x_current-x_p)^2 + (y_current-y_p)^2);
            theta_p = atan2(y_p-y_current, x_p-x_current);
            phi_current = phi_current + U*R*sin(theta_p)/r;
            psi_current = psi_current - U*R*cos(theta_p)/r;
        end
        
        % 更新流线轨迹坐标
        x_streamline(j) = x_current + dx*phi_current/U;
        y_streamline(j) = y_current + dy*psi_current/U;
    end
    
    % 绘制流线
    hold on;
    plot(x_streamline, y_streamline);
end

% 绘制圆柱
hold on;
plot(x, y, 'k');
axis equal;