Demo entry 6358143

aaa

   

Submitted by anonymous on Apr 24, 2017 at 19:35
Language: Matlab. Code size: 5.3 kB.

%Initial Conditions SI unit 
clear; 
t0=0; 
tfinal=0.0004; 
tspan=[0 tfinal]; 
N=100; 
L=0.00319; 
h=0.15; 
%Initial position 
zi=[0 0 0 0 2 0 0 0]; 
for i=1:N 
z1(i,:)=h-(h/N)*i; 
z2(i,:)=0; 
x1(i,:)=10^-3*L*cos(((2*pi)/lambda)*z1(i,:))*((h-z1(i,:))/h); 
x2(i,:)=0; 
zdt1(i,:)=10^-3*L*sin(((2*pi)/lambda)*z1(i,:))*((h-z1(i,:))/h); 
y2(i,:)=0; 
s1(i,:)=0; 
s2(i,:)=0;
z(i,:)=[x1(i,:),x2(i,:),zdt1(i,:),y2(i,:),z1(i,:),z2(i,:),s1(i,:),s2(i,:)]'; 
zi=[zi,z(i,:)]; 
end 
zt0=zi'; 
[t,zt]=ode45('ESP',tspan,zt0); 
%3D plot 
for i=1:8:8*N+8 
A(i)=zt(end,i); 
B(i)=zt(end,i+2); 
C(i)=zt(end,i+4); 
end 
Az=A(1:8:8*N+8); 
Bz=B(1:8:8*N+8); 
Cz=C(1:8:8*N+8); 
plot3(Az,Bz,Cz) 
%Defining function ESP 
function ztdot=project1(t,zt) 
%electrospinning parameters 
m=0.28e-8; 
e=2.83e-9; 
vo=15000; 
h=0.15; 
L=3e-3; 
G=10e4; 
mu=7; 
ao=150e-6; 
alpha=0.7; 
lambda=3e6; 
N=100; %number of beads modelled 
ztadot=[0 0 0 0 0 0 0 0]; 
Rx=zt(1:8:end); 
Ry=zt(3:8:end); 
Rz=zt(5:8:end); 
%Coulombic Force 
for i=1:N+1; 
for j=1:N+1 
if i~=j; 
R(i,j)=sqrt((Rx(i)-Rx(j))^2+(Ry(i)-Ry(j))^2+(Rz(i)-Rz(j))^2); 
Fcx(i,j)=((e^2/R(i,j)^2)*(Rx(i)-Rx(j))); 
Fcy(i,j)=((e^2/R(i,j)^2)*(Ry(i)-Ry(j))); 
Fcz(i,j)=((e^2/R(i,j)^2)*(Rz(i)-Rz(j))); 
end 
end 
end 
for j=1:N+1 
fx(j)=s1m(Fcx(j,:)); 
fy(j)=s1m(Fcy(j,:));
fz(j)=s1m(Fcz(j,:)); 
end 
fc=[fx;zeros(1,N+1);fy;zeros(1,N+1);fz;zeros(3,N+1)]; 
Fc=fc(1:end); 
for i=9:8:8*N+8; 
lu(i)=sqrt((zt(i-8)-zt(i))^2+(zt(i-6)-zt(i+2))^2+(zt(i-4)-zt(i+4))^2); 
au(i)=sqrt((ao^2*L)/lu(i)); 
end 
for i=1:8:8*N; 
ld(i)=sqrt((zt(i)-zt(i+8))^2+(zt(i+2)-zt(i+10))^2+(zt(i+4)-zt(i+12))^2); 
ad(i)=sqrt((ao^2*L)/ld(i)); 
end 
i=1; %first bead 
s1dot(i,:)=0; 
s2dot(i,:)=(G/(ld(i))^2)*(((zt(i)-zt(i+8))*zt(i+1))+((-zt(i)+zt(i+8))*zt(i+9))+((zt(i+2)-zt(i+10))*zt(i+3))+((-zt(i+2)+zt(i+10))*zt(i+11))+((zt(i+4)-zt(i+12))*zt(i+5))+((-zt(i+4)+zt(i+12))*zt(i+13)))-(G/mu)*zt(i+7); 
x1dot(i,:)=zt(i+1); 
x2dot(i,:)=zt(i+2); 
zdt1dot(i,:)=zt(i+3); 
y2dot(i,:)=zt(i+4); 
z1dot(i,:)=zt(i+5); 
z2dot(i,:)=zt(i+6); 
ztdot(1,:)=[x1dot(i,:) x2dot(i,:) zdt1dot(i,:) y2dot(i,:) z1dot(i,:) z2dot(i,:) s1dot(i,:) s2dot(i,:)]; 
for i=9:8:8*N; 
if zt(i+4)<0; 
zt(i+4)=0; 
end 
if zt(i)>0; 
signx=1; 
end 
if zt(i)<0; 
signx=-1; 
end 
if zt(i+2)>0; 
signy=1; 
end 
if zt(i+2)<0; 
signy=-1; 
end 
m1(i)=(zt(i+2)-zt(i+10))/(zt(i)-zt(i+8)); 
m2(i)=(zt(i-6)-zt(i+2))/(zt(i-8)-zt(i)); 
dsxu(i)=abs(zt(i-8)-zt(i)); 
dsyu(i)=abs(zt(i-6)-zt(i+2)); 
dsxd(i)=abs(zt(i)-zt(i+9)); 
dsyd(i)=abs(zt(i+2)-zt(i+10)); 
ds1(i)=sqrt(dsxu(i)^2+dsyu(i)^2); 
ds2(i)=sqrt(dsxd(i)^2+dsyd(i)^2);
k(i)=atan((m2(i)-m1(i))/(1+m1(i)*m2(i)))/(ds1(i)+ds2(i)); 
s1dot(i,:)=(G/(lu(i))^2)*(((zt(i-8)-zt(i))*zt(i-7))+((-zt(i-8)+zt(i))*zt(i+1))+((zt(i-6)-zt(i+2))*zt(i-5))+((-zt(i-6)+zt(i+2))*zt(i+3))+((zt(i-4)-zt(i+4))*zt(i-3))+((-zt(i-4)+zt(i+4))*zt(i+5)))-(G/mu)*zt(i+6); 
s2dot(i,:)=(G/(ld(i))^2)*(((zt(i)-zt(i+9))*zt(i+1))+((-zt(i)+zt(i+8))*zt(i+9))+((zt(i+2)-zt(i+10))*zt(i+3))+((-zt(i+2)+zt(i+10))*zt(i+11))+((zt(i+4)-zt(i+12))*zt(i+5))+((-zt(i+4)+zt(i+12))*zt(i+13)))-(G/mu)*zt(i+7); 
x1dot(i,:)=zt(i+1); 
x2dot(i,:)=(Fc(i)-((pi*(au(i)^2)*zt(i+6))/lu(i))*(zt(i-8)-zt(i))-((pi*(ad(i)^2)*zt(i+7))/ld(i))*(zt(i)-zt(i+8))-((alpha*pi*(((au(i)+ad(i))^2)/4)*k(i))/sqrt(zt(i)^2+zt(i+2)^2))*(abs(zt(i))*signx))/m; 
zdt1dot(i,:)=zt(i+3); 
y2dot(i,:)=(Fc(i+2)-((pi*(au(i)^2)*zt(i+6))/lu(i))*(zt(i-6)-zt(i+2))-((pi*(ad(i)^2)*zt(i+7))/ld(i))*(zt(i+2)-zt(i+10))-((alpha*pi*(((au(i)+ad(i))^2)/4)*k(i))/sqrt(zt(i)^2+zt(i+2)^2))*(abs(zt(i+2))*signy))/m; 
z1dot(i,:)=zt(i+5); 
z2dot(i,:)=(Fc(i+4)-(e*vo)/h+((pi*(au(i)^2)*zt(i+6))/lu(i))*(zt(i-4)-zt(i+4))-((pi*(ad(i)^2)*zt(i+7))/ld(i))*(zt(i+4)-zt(i+12)))/m; 
ztdot(i,:)=[x1dot(i,:) x2dot(i,:) zdt1dot(i,:) y2dot(i,:) z1dot(i,:) z2dot(i,:) s1dot(i,:) s2dot(i,:)]; 
ztadot=[ztadot,ztdot(i,:)]; 
end 
ztmdot=ztadot(:,9:8*N); 
i=8*N+1; %last bead 
if zt(i+4)<0; 
zt(i+4)=0; 
end 
if zt(i)>0; 
signx=1; 
end 
if zt(i)<0; 
signx=-1; 
end 
if zt(i+2)>0; 
signy=1; 
end 
if zt(i+2)<0; 
signy=-1; 
end 
m1(i)=zt(i+2)/zt(i); 
m2(i)=(zt(i-6)-zt(i+2))/(zt(i-8)-zt(i)); 
dsxu(i)=abs(zt(i-8)-zt(i)); 
dsyu(i)=abs(zt(i-6)-zt(i+2)); 
dsxd(i)=abs(zt(i)); 
dsyd(i)=abs(zt(i+2));
ds1(i)=sqrt(dsxu(i)^2+dsyu(i)^2); 
ds2(i)=sqrt(dsxd(i)^2+dsyd(i)^2); 
k(i)=atan((m2(i)-m1(i))/(1+m1(i)*m2(i)))/(ds1(i)+ds2(i)); 
s1dot(i,:)=(G/(lu(i))^2)*(((zt(i-8)-zt(i))*zt(i-7))+((-zt(i-8)+zt(i))*zt(i+1))+((zt(i-6)-zt(i+2))*zt(i-5))+((-zt(i-6)+zt(i+2))*zt(i+3))+((zt(i-4)-zt(i+4))*zt(i-3))+((-zt(i-4)+zt(i+4))*zt(i+5)))-(G/mu)*zt(i+6); 
s2dot(i,:)=0; 
x1dot(i,:)=zt(i+1); 
x2dot(i,:)=(Fc(i+2)-((pi*(au(i)^2)*zt(i+6))/lu(i))*(zt(i-9)-zt(i))-((alpha*pi*(((au(i))^2)/4)*k(i))/sqrt(zt(i)^2+zt(i+2)^2))*(abs(zt(i))*signx))/m; 
zdt1dot(i,:)=zt(i+3); 
y2dot(i,:)=(Fc(i+2)-((pi*(au(i)^2)*zt(i+6))/lu(i))*(zt(i-6)-zt(i+2))-((alpha*pi*(((au(i))^2)/4)*k(i))/sqrt(zt(i)^2+zt(i+2)^2))*(abs(zt(i+2))*signy))/m; 
z1dot(i,:)=zt(i+5); 
z2dot(i,:)=(Fc(i+4)-(e*vo)/h+((pi*(au(i)^2)*zt(i+6))/lu(i))*(zt(i-4)-zt(i+4)))/m; 
ztdot(i,:)=[x1dot(i,:) x2dot(i,:) zdt1dot(i,:) y2dot(i,:) z1dot(i,:) z2dot(i,:) s1dot(i,:) s2dot(i,:)]; 
ztdot=[ztdot(1,:) ztmdot ztdot(8*N+1,:)]';

This snippet took 0.02 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).