Asked by Hollis Williams
on 7 Jul 2019

I have a code with MATLAB where I set the number of boundary nodes around a sphere and then evaluate functions at the nodes to get an expression for drag force which is then divided by the appropriate analytic expression for drag obtained without numerical methods. Obviously this ratio just gives 1.00000 if the number of nodes is sufficiently high, but I would like a graph which shows the convergence as the number of nodes runs from 1 to 100, so I need a code which evaluates through a range of values of N then plots the ratio of numerical result divided by analytic expression on a graph.

Let me know if this is not clear or if I need to post some of the code. The code is a bit of a mess, so I will just explain the key. N is the number of nodes, each node is parametrised by the angle phi and the angle theta (which I call alpha so as not to confuse with temperature). The number of nodes is always N_phi*N_alpha, where N_phi = N_alpha, so for example I need to run through the case where they are both 2 (ie. so that N=4), both 3, etc and so up to where they are both 10, and then plot these cases on a graph of number of nodes on the sphere vs accuracy of the result (this is just the ratio of the numerical result divided by the analytic result, which is ratio1 in the code). The result is calculated with a numerical singularity method where you have singularities and boundary nodes and apply functions which take the radial distance between a singularity and a boundary node for every node on the boundary to create a linear system corresponding to the point forcing on every node, which is then solved for an overall drag force.

clear all;

close all;

clc;

VZ=1;

r=1;

Kn=0.7*sqrt(2/pi);

rsin=0.2;

H=ones(3,3)-eye(3);

int=0;

Int=0;

INT=0;

N_alpha=10;

N_phi=10;

N = (N_phi .* N_alpha);

M(3*N,3*N)=0;

for alphaindex=1:N_alpha

alpha=(alphaindex-1)*2*pi/N_alpha; %index positions of nodes in spherical polars%

for phiindex=1:N_phi

phi=(phiindex)*pi/(N_phi+1);

int=int+1;

x(alphaindex,phiindex)=r*cos(alpha)*sin(phi);

y(alphaindex,phiindex)=r*sin(alpha)*sin(phi); %convert to Cartesian coordinates%

z(alphaindex,phiindex)=r*cos(phi);

b(1:3,int)=[x(alphaindex,phiindex) y(alphaindex,phiindex) z(alphaindex,phiindex)];

end

end

for alphaindex=1:N_alpha

alpha=(alphaindex-1)*2*pi/N_alpha; %repeat for the singularities%

for phiindex=1:N_phi

phi=(phiindex)*pi/(N_phi+1);

Int=Int+1;

sx(alphaindex,phiindex)=rsin*cos(alpha)*sin(phi);

sy(alphaindex,phiindex)=rsin*sin(alpha)*sin(phi);

sz(alphaindex,phiindex)=rsin*cos(phi);

sing(1:3,Int)=[sx(alphaindex,phiindex) sy(alphaindex,phiindex) sz(alphaindex,phiindex)];

end

end

A = b./sqrt(sum(b.*b));

B = cross(A,repmat([0 0 1]',1,size(A,2))); %obtain normals and tangents at all nodes%

C = cross(A,B);

D = [0 0 1]*A;

E = [0 0 1]*B; %contract with velocity on RHS in Cartesian%

F = [0 0 1]*C;

V = [D;E;F];

VV = reshape(V,[],1); %new velocity at nodes in normal-tangentials%

for Ind=1:N

for IndS=1:N

unit=[0 0 1]';

normal=(b(1:3,Ind)/(norm(b(1:3,Ind))))';

t = (cross((b(1:3,Ind)/(norm(b(1:3,Ind)))),unit))';

s = cross((b(1:3,Ind)/(norm(b(1:3,Ind)))),(cross((b(1:3,Ind)/(norm(b(1:3,Ind)))),unit)))';

rij=(b(1:3,Ind)-sing(1:3,IndS))';

M(1+3*(Ind-1),1+3*(IndS-1))=PartA(rij,normal);

M(1+3*(Ind-1),2+(3*(IndS-1)))=PartB(rij,normal); %create matrix for linear system%

M(1+3*(Ind-1),3+(3*(IndS-1)))=PartC(rij,normal);

M(2+3*(Ind-1),1+3*(IndS-1))=PartD(rij,t);

M(2+3*(Ind-1),2+3*(IndS-1))=PartE(rij,t);

M(2+3*(Ind-1),3+3*(IndS-1))=PartF(rij,t);

M(3+3*(Ind-1),1+3*(IndS-1))=PartG(rij,s);

M(3+3*(Ind-1),2+3*(IndS-1))=PartH(rij,s);

M(3+3*(Ind-1),3+3*(IndS-1))=PartI(rij,s);

end

end

f = M\VV; %solve linear system%

dragFX=0;dragFY=0; dragFZ=0; %calculate drag force%

for ind=1:N

dragFX=f(1+(ind-1)*3)+dragFX;

dragFY=f(2+(ind-1)*3)+dragFY;

dragFZ=f(3+(ind-1)*3)+dragFZ;

end

dragFX=-1*dragFX; dragFY=-1*dragFY; dragFZ=-1*dragFZ;

stokesDrag=-6*pi*VZ*Kn;

slipDrag=stokesDrag*(1+2*Kn/sqrt(2/pi))/(1+3*Kn/sqrt(2/pi)); %compare with analytic solution%

ratio1 = dragFZ/stokesDrag;

ratio2=dragFZ/slipDrag;

Opportunities for recent engineering grads.

Apply Today
## 4 Comments

## darova (view profile)

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/470582-graph-of-convergence-for-drag-force#comment_722075

## Hollis Williams (view profile)

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/470582-graph-of-convergence-for-drag-force#comment_722124

## darova (view profile)

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/470582-graph-of-convergence-for-drag-force#comment_722130

## Hollis Williams (view profile)

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/470582-graph-of-convergence-for-drag-force#comment_722139

Sign in to comment.