Graph of Convergence for Drag Force

2 views (last 30 days)
Hollis Williams
Hollis Williams on 7 Jul 2019
Commented: 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;
  4 Comments
Hollis Williams
Hollis Williams on 7 Jul 2019
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.

Sign in to comment.

Answers (0)

Tags

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!