How do determine the pressure? Please check me my code.

if true
  % code
end
t = 0:0.1:10*pi;
r = linspace (0, 1, numel (t));
z = linspace (0, 1, numel (t));
a = 0.009;
Gamm_inf = 2000;
nu = 3;
p0 = 1;
ro = 2;
vth = (Gamm_inf)./(2*pi*nu.*r);
p_r = p0 + ro.*int(vth*vth./r,r,0,r)-(ro*a^2)/2.*(r^2 + 4.*z^2);
plot(p_r,r);
plot(p_r,z);
plot3(p_r,r,z);

Answer 1

[Star Strider]

To do numerical integration, you need to use the integral function. Here, it is also necessary to use arrayfun to integrate over your ‘r’ vector. I changed ‘vth’ and also the plots, so ‘p_r’ is plotted as a function of ‘r’ (and ‘z’).

This runs, without errors or warnings. I will let your determine if it is correct:

t = 0:0.1:10*pi;
r = linspace (0, 1, numel (t));
z = linspace (0, 1, numel (t));
a = 0.009;
Gamm_inf = 2000;
nu = 3;
p0 = 1;
ro = 2;
vth = @(r) (Gamm_inf)./(2*pi*nu.*r) .* (1 - exp(-(a*r/(2*nu))));
p_r = p0 + ro.*arrayfun(@(r)integral(@(r)vth(r).*vth(r)./r,1E-8,r, 'ArrayValued',1),r)-(ro*a^2)/2.*(r.^2 + 4.*z.^2);

figure
plot(r, p_r);
grid

figure
plot3(r, z ,p_r);
grid on

Experiment to get the result you want.