How can I appear the streamline around the spheres?

Question

Shreen El-Sapa on 22 Apr 2024

0
Link

Direct link to this question

https://in.mathworks.com/matlabcentral/answers/2110391-how-can-i-appear-the-streamline-around-the-spheres

Edited: Voss on 22 Apr 2024

stream1.m

Open in MATLAB Online

clc

A=[ -0.2

0.1

-0.4

1.0

-2.1

3.9

-6.9

11.6

-18.7

29.2

-44.0

64.9

-93.5

132.6

-184.2

253.6

-343.1

462.3

-613.4

815.3

-1068.2

1413.3

-1842.6

2462.0

-3232.1

4517.3

-6127.8

10894.9

-17024.9

9869.9];

B=[ 279.2

34.9

-96.3

177.4

-248.7

283.6

-271.0

223.3

-160.5

102.6

-58.7

30.5

-14.4

6.3

-2.5

0.9

-0.3

0.1

0.0

0.0];

C=[ -4.6

-0.3

0.3

-0.2

0.1

0.0

0.0];

AA=[ -1.4

-0.1

0.4

-1.1

3.0

-7.5

17.7

-40.2

88.1

-188.2

390.9

-797.9

1594.7

-3150.4

6115.4

-11791.2

22393.6

-42435.9

79327.3

-148756.4

275335.2

-515215.3

950981.9

-1800517.1

3351792.8

-6647748.5

12804165.0

-32341874.0

71833152.0

-59216272.0];

BB=[ 15518.5

-406.6

1241.4

-2229.4

3366.8

-4306.2

4732.0

-4611.8

4004.3

-3152.0

2255.2

-1486.4

902.0

-510.4

268.4

-132.9

61.6

-27.2

11.3

-4.5

1.7

-0.6

0.2

-0.1

0.0

0.0];

CC=[ -15.6

0.3

-0.4

0.3

-0.2

0.1

0.0

0.0];

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

a = 1 ; %RADIUS

L=.22;

etta=0.01; %0.02;0.05;0.1d0

u2=1; delta=1.5; ettap=.01; alpha=2.d0; %C=a1/a2=0.1

xi=1./sqrt(etta);xi1=1./sqrt(ettap);

alpha1=sqrt((xi.^2+sqrt(xi.^4-4.*xi.^2.*alpha.^2))./sqrt(2));

alpha2=sqrt((xi.^2-sqrt(xi.^4-4.*xi.^2.*alpha.^2))./sqrt(2));

dd=5;

c =-a/L;

b =a/L;

m =a*40; % NUMBER OF INTERVALS

[x,y]=meshgrid([c+dd:(b-c)/m:b],[c:(b-c)/m:b]');

[I J]=find(sqrt(x.^2+y.^2)<(a-.0));

if ~isempty(I);

x(I,J) = 0; y(I,J) = 0;

end

r=sqrt(x.^2+y.^2);

t=atan2(y,x);

r2=sqrt(r.^2+dd.^2-2.*r.*dd.*cos(t));

zet=(r.^2-r2.^2-dd.^2)./(2.*r2.*dd);

%for i1=1:length(x);

% for k1=1:length(x);

% if sqrt(x(i1,k1).^2+y(i1,k1).^2)>1./L;

% r(i1,k1)=0;r2(i1,k1)=0;

% end

%end

warning off

psi1=0;

for i=2:7

Ai=A(i-1);Bi=B(i-1);Ci=C(i-1);AAi=AA(i-1);BBi=BB(i-1);CCi=CC(i-1);

psi1=psi1+(Ai.*r.^(-i+1)+r.^(1./2).*besselk(i-1./2,r.*alpha1).*Bi+r.^(1./2).*besselk(i-1./2,r.*alpha2).*Ci).*gegenbauerC(i,-1./2, cos(t))+(AAi.*r2.^(-i+1)+r2.^(1./2).*besselk(i-1./2,r2.*alpha1).*BBi+r2.^(1./2).*besselk(i-1./2,r2.*alpha2).*CCi).*gegenbauerC(i,-1./2,zet);

end

hold on

%[DH1,h1]=contour(x,y,psi1,25,'-k','LineWidth',1.1); %,psi2,'--k',psi2,':k'

%[DH1,h1]=contour(x,y,psi1);

%p1=contour(x,y,psi1,[0.3 0.3],'k','LineWidth',1.1); %,'ShowText','on'

%p2=contour(x,y,psi1,[0.4 0.4],'r','LineWidth',1.1);

%p3=contour(x,y,psi1,[0.5 0.5],'g','LineWidth',1.1);

%p4=contour(x,y,psi1,[0.6 0.6],'b','LineWidth',1.1);

%p5=contour(x,y,psi1,[0.7 0.7],'c','LineWidth',1.1);

%p6=contour(x,y,psi1,[0.8 0.8],'m','LineWidth',1.1);

%p7=contour(x,y,psi1,[0.9 0.9],'y','LineWidth',1.1);

p1=contour(x,y,psi1,[-0.001 0.001],'k','LineWidth',1.1); %,'ShowText','on'

p2=contour(x,y,psi1,[-0.005 .005],'r','LineWidth',1.1);

%p3=contour(x,y,psi1,[0.1 0.1],'g','LineWidth',1.1);

%p4=contour(x,y,psi1,[0.4 0.4],'b','LineWidth',1.1);

%p5=contour(x,y,psi1,[0.6 0.6],'c','LineWidth',1.1);

%p6=contour(x,y,psi1,[0.8 0.8],'m','LineWidth',1.1);

%%%%%%%%%%%%%%% $\frac{\textstyle a_1+a_2}{\textstyle h}=6.0,\;

hold on

t3 = linspace(0,2*pi,1000);

h2=0;

k2=0;

rr2=2;

x2 = rr2*cos(t3)+h2;

y2 = rr2*sin(t3)+k2;

set(plot(x2,y2,'-k'),'LineWidth',1.1);

fill(x2,y2,'w')

hold on

t2 = linspace(0,2*pi,1000);

h=dd;

k=0;

rr=1;

x1 = rr*cos(t2)+h;

y1 = rr*sin(t2)+k;

set(plot(x1,y1,'-k'),'LineWidth',1.1);

fill(x1,y1,'w')

%axis square;

axis('equal')

axis on

%xticklabels([])

%yticklabels([])

%legend('0.01','0.05','0.1','0.4','0.6','0.8','Location','northwest')

%title('$\frac{\beta_1}{a_1\mu}=\frac{a_1\beta_2}{\mu}=1.0,\;R_{H}=1.0,\;\frac{a_2}{a_1}=2.0$','Interpreter','latex','FontSize',12,'FontName','Times New Roman','FontWeight','Normal')

%title('$(a)\;\; R_{H}=1.0,\;\frac{\kappa}{\mu}=4.0$','Interpreter','latex','FontSize',12,'FontName','Times New Roman','FontWeight','Normal')

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

view([90 90])