Graph of Convergence for Drag Force

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;

[darova]

POst your code please. Maybe you have some scheme or picture how it looks like?

[Hollis Williams]

Hi darova, I've posted a code. It's a bit of a mess, but the picture is basically that you have a velcoity at every node on a surface and then singular solutions placed at singularity sites inside the surface and that you then work with the distance between the nodes and the singularity sites to get a point forcing at every boundary node (basically just application of Green's function) and then solve the linear system to get a drag force result.

So essentially, I just want to cycle through the possible values which the number of nodes can take and then plot the accurary of the result to show that the result converges to the analytic result as the number of nodes is increased to 100. The complication is that the number of nodes is defined by two separate integer parameters N_alpha and N_phi and I essentially need to run through the result for all the cases where N_alpha and N_phi are equal, so starting with N_alpha = N_phi =1, as far as N_alpha=N_phi=10, which gives 100 nodes.

Let me know if you require more explanation, as this stuff seems to be a bit difficult to explain to people unfortunately and I am trying to improve my ability to explain it.

[darova]

Ok we have sphere, nodes, normals. Am i right?

[Hollis Williams]

That is correct. We have a flow on the sphere then we split the components of the force into the normal and tangential components, although this is not necessary.

You index the nodes on the sphere in spherical polar coordinates and the singularity sites, which have the same positions but different radial distances as they are on a smaller sphere inside the sphere. The positions of the nodes and singularites are converted to Cartesian coordinates. We obtain a vector which gives the velocity due to the flow on the surface at each node in normal-tangential coordinates (not necessary to have these coordinates, but an exercise to try a coordinate conversion).

You have functions corresponding to particular singular solutions (the solution is one called the Stokeslet for this case). You apply the singular solution for every boundary node but it's singular so what you do is loop twice in the code and then evaluate that function for each node for the distance from the node to each singularity, so for the first node you add up the function evaluating for the distance from that node to the first singularity plus evaluation for the distance from the same node to the next singularity.

Then you do that for the next node along. This all gives you a linear system where the unknown is a force vector f and the right-side is the velocity we mentioned for the flow on the wall. You then use that to calculate a total drag force on the sphere, which is compared to an analytic result. I am then essentially just running the process multiple times at once but running through the possible cases where N_alpha and N_phi have the same values, then plot those cases for accuracy vs N to show convergence.

If this doesn't make sense, perhaps I could link to the article with the calculations I am following and hopefully that makes more sense.