Warning: Integer operands are required for colon operator when used as index.
> In Case1_Ijuptilk_130822 (line 75)
Warning: Integer operands are required for colon operator when used as index.
> In Case1_Ijuptilk_130822 (line 83)
%% initialization
clear
clc
close all
%% Model parameters
global k1 eta1 alpha3 gamma1 d N Phi xi h A B C G
k1 = 6*10^(-12); % elastic constant
eta1 = 0.0240; % viscosity
xi = 0.585; % activity parameter
alpha3 = -0.001104; % viscosity
gamma1 = 0.1093; % viscosity
Theta = 0.0001;
d = 0.0002;
% Plotting r(q)
alpha=1-alpha3^2/gamma1/eta1;
q=0:0.01:(5*pi);
r=q-(1-alpha)*tan(q)+(alpha3*xi*alpha/eta1)./(4*k1*q.^2/d^2-alpha3*xi/eta1).*tan(q);
figure
plot(q,r)
axis([0 5*pi -20 20])
% Simplify model parameters
A = k1/gamma1;
B = alpha3/gamma1;
C = eta1 - alpha3*B;
G = alpha3*A;
N = 30; % N must be even
%% Numerical setup
% step size
h = d/N;
% t span
tspan = [0:0.1:100];
nt=length(tspan);
% range of z
%z=linspace(0,d/2,N+1)
z=linspace(0,d,N+1);
% initial conditions
theta0 = Theta*sin(pi*z/d);
v0 = zeros(1,N+1);
theta0_int=theta0(2:N);
v0_int=v0(2:N);
u0 = [theta0_int'; v0_int'];
% Constant mass matrix M: We use M to represent the left hand side of the system of equations.
M1=eye(N-1,N-1);
M2=zeros(N-1,N-1);
M=[M1 M2;M2 M2];
% Boundary Conditions
Phi = 0;
%% ode solver
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',1e-6);
[t,y] = ode15s(@lcode1,tspan,u0, options);
% Extract the solution for theta and v
theta = [Phi*ones(length(t), 1) y(:,1:N-1) Phi*ones(length(t), 1)];
v = [zeros(length(t), 1) y(:,N:(2*N-2)) zeros(length(t), 1)];
% theta at the middle of the layer (i.e., z=d/2)
theta_middle = theta(:, N/2);
%% Plot the solution.
figure
subplot(2,1,1)
plot(z, theta(1:(nt/10):nt,:))
xlabel('z(\mum)')
ylabel('\theta(z,t)(rad)')
legend('t_0 =0', 't_1=10', 't_2=20', 't_3=30', 't_4=40', 'Location','bestoutside')
%title('Director angle','FontSize',10)
%[xmin(z) xmax(z) ymin(theta) ymax(theta)];
subplot(2,1,2)
plot(z, v(1:(nt/10):nt,:))
xlabel('z(\mum)')
ylabel('v(z,t)(m/s)')
legend('t_0 =0', 't_1=10', 't_2=20', 't_3=30', 't_4=40', 'Location','bestoutside')
%hF=gcf;
%exportgraphics(hF,'thetav007neg.png','Resolution',300)
%title('Flow velocity')
%[xmin(z) xmax(z) ymin(theta) ymax(theta)];
figure
plot(t,theta_middle)
xlabel('t (s)')
ylabel('\theta(d/2,t)(rad)')
%title('Director angle at the middle of the layer', 'FontSize',6)
%hF=gcf;
%exportgraphics(hF,'thetamiddleneg.png','Resolution',300)
%% Functions
% Right hand side of the ODEs: F(U)
function rhsode = lcode1(t,y)
global Phi xi h A B C G N
% initialize theta and v
theta = y(1:(N-1));
v = y(N:(2*N-2));
% theta equations
% theta equations
% for the positon to the right of the left hand boundary z=0
rhsode(1,1) = (A/(h^2))*(theta(2)-2*theta(1)+ Phi) - (B/(2*h))*(v(2)-0);
% for all other internal positions 0<z<d
for i=2:(N-2)
rhsode(i,1) = (A/(h^2))*(theta(i+1)-2*theta(i)+theta(i-1)) - (B/(2*h))*(v(i+1)-v(i-1));
end
% for the positon to the left of the right hand boundary z=d
rhsode(N-1,1) = (A/(h^2))*(Phi-2*theta(N-1)+ theta(N-2)) - (B/(2*h))*(0-v(N-2));
% v equations [REMEMBER theta the RHS index for the v equations is N-1 more than the indices of the theta and v variables in this function]
% for the two positons to the right of the left hand boundary z=0
rhsode(N,1) = (G/(h^3))*(-Phi + 3*theta(1) - 3*theta(2) + theta(3)) + (C/(h^2))*(v(2) -2*v(1)+ 0) + (xi/(2*h))*(theta(2)-Phi);
rhsode(N+1,1) = (G/(2*h^3))*(Phi - 2*theta(1) + 2*theta(3)- theta(4)) + (C/(h^2))*(v(3) -2*v(2) + v(1)) + (xi/(2*h))*(theta(3)-theta(1));
% for all other internal positions 0<z<d
for i=(N+2):(2*N-4)
rhsode(i,1) = (G/(2*h^3))*(theta(i-2-(N-1)) - 2*theta(i-1-(N-1)) + 2*theta(i+1-(N-1))- theta(i+2-(N-1))) +(C/(h^2))*(v(i+1-(N-1)) -2*v(i-(N-1)) + v(i-1-(N-1))) + (xi/(2*h))*(theta(i+1-(N-1))-theta(i-1-(N-1)));
end
% for the two positons to the left of the right hand boundary z=d
rhsode(2*N-3,1) = (G/(2*h^3))*(theta(N-4) - 2*theta(N-3) + 2*theta(N-1) - Phi) +(C/(h^2))*(v(N-1) -2*v(N-2) + v(N-3)) + (xi/(2*h))*(theta(N-1)-theta(N-3));
rhsode(2*N-2,1) = (G/(h^3))*(-theta(N-3) + 3*theta(N-2) - 3*theta(N-1) + Phi) +(C/(h^2))*(0 -2*v(N-1) + v(N-2)) + (xi/(2*h))*(Phi-theta(N-2));
end
