Graph of Convergence for Drag Force
1 visualización (últimos 30 días)
Mostrar comentarios más antiguos
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 comentarios
Respuestas (0)
Ver también
Categorías
Más información sobre Calendar en Help Center y File Exchange.
Productos
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!