# Demo entry 6358144

aaa

Submitted by aa on Apr 24, 2017 at 19:38
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);
end
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);
zdt1dot(i,:)=zt(i+3);
z1dot(i,:)=zt(i+5);
ztdot(i,:)=[x1dot(i,:) x2dot(i,:) zdt1dot(i,:) y2dot(i,:) z1dot(i,:) z2dot(i,:) s1dot(i,:) s2dot(i,:)];
end
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.