How to solve system of PDE’s using crank nicolson method to get graphical interpretations of equations? How to get skin friction and nusselt number using this code?
i have attached my matlb code for solving system of PDE ∂u/∂t=(∂^2 u)/〖∂η〗^2 +α (∂^3 u)/〖∂η〗^3 +δcφ+Grθ-Mu
∂θ/∂t=1/Pr (∂^2 θ)/〖∂η〗^2 +D_f (∂^2 φ)/〖∂η〗^2
∂φ/∂t=1/S_c (∂^2 φ)/〖∂η〗^2 +S_r (∂^2 θ)/〖∂η〗^2
when η→0
when η→0:u=sin(wt),θ=1 ,φ=1
when η→∞:u→0,θ→0,φ→0
But,i’m confused how to set boundary conditions in this code. Also i want to set these skin friction and nusselt number expressions: C_f=-(∂u/∂η) at η=0
Nu〖Re〗_x^(-1)=-(∂θ/∂η) at η=0
in this code;
function yy= crank
n=105;
h=0.01;
m=0.1;
B=0; %B=0.0001;
Gr=2; %Gr=0.1=Delt. Gs=Dels
A=2; %A=alpha
Pr=5;
Gs=1; %3.5
Sc=1.5;
M=9; %M=30,12
Df=5;
Sr=10;
Ec=0.1; % 0.01-0.1(0-1)
t(1)=0;
s(1)=0;
for j=2:n
t(j)=t(j-1)+h;
end
for j=2:n
s(j)=s(j-1)+m;
end
k=1;
for j=1:n
aa(j,k)=exp(-0.15*Pr^2*t(j)*Ec^-0.1/Gr^0.1); % Guess M*Gr*Gs*S
bb(j,k)=exp(-t(j)*Sc/(Ec)^0.01*M); %bb(j,k)=exp(-Pr^2*t(j)/Ec);
cc(j,k)=exp(-Sc*2*t(j)/Ec^0.3);
end
a{1,k}=[1 0 0;0 1 0;0 0 1];
for j=2:n-1
a{j,k}=[-(1/h)-(1/h^2)-2*A-M Gr Gs;0 -(1/h)-(1/Pr*(h^2)) -Df/h^2;0 -Sr/h^2 -(1/h)-(1/Sc*h^2)];
end
a{n,k}=[1 0 0;0 1 0;0 0 1];
for j=2:n-1
b{j,k}=[(1/2*h^2)+(A/h^3) 0 0;0 1/2*Pr*h^2 Df/2*h^2;0 Sr/2*h^2 1/2*Sc*h^2]; % Lower diagonal entries
end
b{n,k}=[0 0 0;0 0 0;0 0 0];
c{1,k}=[0 0 0;0 0 0;0 0 0];
for j=2:n-1
c{j,k}=[(1/2*h^2)+(A/h^3) 0 0;0 1/2*Pr*h^2 Df/2*h^2;0 Sr/2*h^2 1/2*Sc*h^2];
end
r1(1,k)=0;
r2(1,k)=0;
r3(1,k)=0;
for j=2:n-1
r1(j,k)=-(exp(t(j))/h)-(1/2*h^2)*(1)*(exp(-t(j+1))-2*exp(-t(j))+exp(-t(j-1)))-Gr*exp(-t(j))-Gs*exp(t(j))+M*exp(-t(j))-(A/2*h^3)*(-exp(-t(j+1))+2*exp(-t(j))-exp(-t(j-1)))…
+(aa(j,k)/h)-(1/2*h^2)*(1)*(aa(j+1,k)-2*aa(j,k)+aa(j-1,k))-(A/2*h^3)*( aa(j+1,k)^2-2*aa(j+1,k)*aa(j,k))-Gr*bb(j,k)-Gs*cc(j,k)+M*aa(j,k);
r2(j,k)=(exp((-t(j)))/h)-(1/2*Pr*h^2)*(exp(-t(j+1))-2* exp(-t(j))+exp(-t(j-1)))-Df*(1/2*h^2)*(exp(-t(j+1))-2* exp(-t(j))…
+exp(-t(j-1)))+(bb(j,k)/h)-(1/2*Pr*h^2)*(bb(j+1,k)-2*bb(j,k)+bb(j-1,k))-(Df/2*h^2)*(cc(j+1,k)-2*cc(j,k)+cc(j-1,k));
r3(j,k)=-(exp((-t(j)))/h)-(1/2*Sc*h^2)*(exp(-t(j+1))-2*exp(-t(j))+exp(-t(j-1)))-(Sr/2*h^2)*(exp(-t(j+1))-2*exp(-t(j))+exp(-t(j-1)))+ (cc(j,k)/h)…
-(1/2*Sc*h^2)*(cc(j+1,k)-2*cc(j,k)+cc(j-1,k))-(Sr/2*h^2)*(bb(j+1,k)-2*bb(j,k)+bb(j-1,k));
end
r1(n,k)= 0;
r2(n,k)=0;
r3(n,k)=0;
for j=1:n
rr{j,k}=[r1(j,k);r2(j,k);r3(j,k)];
end
gamma{1,k}=inv(a{1,k})*c{1,k};
for j=2:n-1
a{j,k}=a{j,k}-(b{j,k}*gamma{j-1,k});
gamma{j,k}=inv(a{j,k})*c{j,k};
end
y{1}=inv(a{1})*rr{1};
for j=2:n
y{j}=inv(a{j})*(rr{j}-b{j}*y{j-1});
end
x{n}=y{n};
for j=n-1:-1:1
x{j}=y{j}-(gamma{j})*x{j+1};
end
for j=n:-1:1
u(j,k)=x{j}(1,1);
z(j,k)=x{j}(2,1);
q(j,k)=x{j}(3,1);
end
for j=1:n
xx(j,k)= aa(j,k)+u(j,k);
yy(j,k)= bb(j,k)+z(j,k);
zz(j,k)= cc(j,k)+q(j,k);
end
for j=1:n-1
%U(j)=-((1+(1/B))*((xx(j+1)-xx(j))/h))
%U(j)=-((yy(j+1)-yy(j))/2*h)
end
%ee=meshgrid(xx);
%V=trapz(ee);
%eee=meshgrid(V);
%[XX,YY]=meshgrid(t,s);
%XXX=meshgrid(xx);
%surf(XX,YY,XXX)
%contour(eee)
figure(1)
plot(t,xx,’LineWidth’,2,’MarkerSize’,16,’linestyle’,’-‘,’color’,’k’)
xlabel(‘eta’,’Interpreter’,’tex’,’FontSize’,16,’FontWeight’,’bold’);
ylabel(‘u’,’Interpreter’,’tex’,’FontSize’,16,’FontWeight’,’bold’);
grid on
hold on
figure(2)
plot(t,yy,’LineWidth’,2,’MarkerSize’,16,’linestyle’,’-‘,’color’,’k’)
xlabel(‘eta’,’Interpreter’,’tex’,’FontSize’,16,’FontWeight’,’bold’);
ylabel(‘Theta(eta)’,’Interpreter’,’tex’,’FontSize’,16,’FontWeight’,’bold’);
grid on
hold on
figure(3)
plot(t,zz,’LineWidth’,2,’MarkerSize’,16,’linestyle’,’–‘,’color’,’r’)
xlabel(‘eta’,’Interpreter’,’tex’,’FontSize’,16,’FontWeight’,’bold’);
ylabel(‘Phi(eta)’,’Interpreter’,’tex’,’FontSize’,16,’FontWeight’,’bold’);
grid on
hold oni have attached my matlb code for solving system of PDE ∂u/∂t=(∂^2 u)/〖∂η〗^2 +α (∂^3 u)/〖∂η〗^3 +δcφ+Grθ-Mu
∂θ/∂t=1/Pr (∂^2 θ)/〖∂η〗^2 +D_f (∂^2 φ)/〖∂η〗^2
∂φ/∂t=1/S_c (∂^2 φ)/〖∂η〗^2 +S_r (∂^2 θ)/〖∂η〗^2
when η→0
when η→0:u=sin(wt),θ=1 ,φ=1
when η→∞:u→0,θ→0,φ→0
But,i’m confused how to set boundary conditions in this code. Also i want to set these skin friction and nusselt number expressions: C_f=-(∂u/∂η) at η=0
Nu〖Re〗_x^(-1)=-(∂θ/∂η) at η=0
in this code;
function yy= crank
n=105;
h=0.01;
m=0.1;
B=0; %B=0.0001;
Gr=2; %Gr=0.1=Delt. Gs=Dels
A=2; %A=alpha
Pr=5;
Gs=1; %3.5
Sc=1.5;
M=9; %M=30,12
Df=5;
Sr=10;
Ec=0.1; % 0.01-0.1(0-1)
t(1)=0;
s(1)=0;
for j=2:n
t(j)=t(j-1)+h;
end
for j=2:n
s(j)=s(j-1)+m;
end
k=1;
for j=1:n
aa(j,k)=exp(-0.15*Pr^2*t(j)*Ec^-0.1/Gr^0.1); % Guess M*Gr*Gs*S
bb(j,k)=exp(-t(j)*Sc/(Ec)^0.01*M); %bb(j,k)=exp(-Pr^2*t(j)/Ec);
cc(j,k)=exp(-Sc*2*t(j)/Ec^0.3);
end
a{1,k}=[1 0 0;0 1 0;0 0 1];
for j=2:n-1
a{j,k}=[-(1/h)-(1/h^2)-2*A-M Gr Gs;0 -(1/h)-(1/Pr*(h^2)) -Df/h^2;0 -Sr/h^2 -(1/h)-(1/Sc*h^2)];
end
a{n,k}=[1 0 0;0 1 0;0 0 1];
for j=2:n-1
b{j,k}=[(1/2*h^2)+(A/h^3) 0 0;0 1/2*Pr*h^2 Df/2*h^2;0 Sr/2*h^2 1/2*Sc*h^2]; % Lower diagonal entries
end
b{n,k}=[0 0 0;0 0 0;0 0 0];
c{1,k}=[0 0 0;0 0 0;0 0 0];
for j=2:n-1
c{j,k}=[(1/2*h^2)+(A/h^3) 0 0;0 1/2*Pr*h^2 Df/2*h^2;0 Sr/2*h^2 1/2*Sc*h^2];
end
r1(1,k)=0;
r2(1,k)=0;
r3(1,k)=0;
for j=2:n-1
r1(j,k)=-(exp(t(j))/h)-(1/2*h^2)*(1)*(exp(-t(j+1))-2*exp(-t(j))+exp(-t(j-1)))-Gr*exp(-t(j))-Gs*exp(t(j))+M*exp(-t(j))-(A/2*h^3)*(-exp(-t(j+1))+2*exp(-t(j))-exp(-t(j-1)))…
+(aa(j,k)/h)-(1/2*h^2)*(1)*(aa(j+1,k)-2*aa(j,k)+aa(j-1,k))-(A/2*h^3)*( aa(j+1,k)^2-2*aa(j+1,k)*aa(j,k))-Gr*bb(j,k)-Gs*cc(j,k)+M*aa(j,k);
r2(j,k)=(exp((-t(j)))/h)-(1/2*Pr*h^2)*(exp(-t(j+1))-2* exp(-t(j))+exp(-t(j-1)))-Df*(1/2*h^2)*(exp(-t(j+1))-2* exp(-t(j))…
+exp(-t(j-1)))+(bb(j,k)/h)-(1/2*Pr*h^2)*(bb(j+1,k)-2*bb(j,k)+bb(j-1,k))-(Df/2*h^2)*(cc(j+1,k)-2*cc(j,k)+cc(j-1,k));
r3(j,k)=-(exp((-t(j)))/h)-(1/2*Sc*h^2)*(exp(-t(j+1))-2*exp(-t(j))+exp(-t(j-1)))-(Sr/2*h^2)*(exp(-t(j+1))-2*exp(-t(j))+exp(-t(j-1)))+ (cc(j,k)/h)…
-(1/2*Sc*h^2)*(cc(j+1,k)-2*cc(j,k)+cc(j-1,k))-(Sr/2*h^2)*(bb(j+1,k)-2*bb(j,k)+bb(j-1,k));
end
r1(n,k)= 0;
r2(n,k)=0;
r3(n,k)=0;
for j=1:n
rr{j,k}=[r1(j,k);r2(j,k);r3(j,k)];
end
gamma{1,k}=inv(a{1,k})*c{1,k};
for j=2:n-1
a{j,k}=a{j,k}-(b{j,k}*gamma{j-1,k});
gamma{j,k}=inv(a{j,k})*c{j,k};
end
y{1}=inv(a{1})*rr{1};
for j=2:n
y{j}=inv(a{j})*(rr{j}-b{j}*y{j-1});
end
x{n}=y{n};
for j=n-1:-1:1
x{j}=y{j}-(gamma{j})*x{j+1};
end
for j=n:-1:1
u(j,k)=x{j}(1,1);
z(j,k)=x{j}(2,1);
q(j,k)=x{j}(3,1);
end
for j=1:n
xx(j,k)= aa(j,k)+u(j,k);
yy(j,k)= bb(j,k)+z(j,k);
zz(j,k)= cc(j,k)+q(j,k);
end
for j=1:n-1
%U(j)=-((1+(1/B))*((xx(j+1)-xx(j))/h))
%U(j)=-((yy(j+1)-yy(j))/2*h)
end
%ee=meshgrid(xx);
%V=trapz(ee);
%eee=meshgrid(V);
%[XX,YY]=meshgrid(t,s);
%XXX=meshgrid(xx);
%surf(XX,YY,XXX)
%contour(eee)
figure(1)
plot(t,xx,’LineWidth’,2,’MarkerSize’,16,’linestyle’,’-‘,’color’,’k’)
xlabel(‘eta’,’Interpreter’,’tex’,’FontSize’,16,’FontWeight’,’bold’);
ylabel(‘u’,’Interpreter’,’tex’,’FontSize’,16,’FontWeight’,’bold’);
grid on
hold on
figure(2)
plot(t,yy,’LineWidth’,2,’MarkerSize’,16,’linestyle’,’-‘,’color’,’k’)
xlabel(‘eta’,’Interpreter’,’tex’,’FontSize’,16,’FontWeight’,’bold’);
ylabel(‘Theta(eta)’,’Interpreter’,’tex’,’FontSize’,16,’FontWeight’,’bold’);
grid on
hold on
figure(3)
plot(t,zz,’LineWidth’,2,’MarkerSize’,16,’linestyle’,’–‘,’color’,’r’)
xlabel(‘eta’,’Interpreter’,’tex’,’FontSize’,16,’FontWeight’,’bold’);
ylabel(‘Phi(eta)’,’Interpreter’,’tex’,’FontSize’,16,’FontWeight’,’bold’);
grid on
hold on i have attached my matlb code for solving system of PDE ∂u/∂t=(∂^2 u)/〖∂η〗^2 +α (∂^3 u)/〖∂η〗^3 +δcφ+Grθ-Mu
∂θ/∂t=1/Pr (∂^2 θ)/〖∂η〗^2 +D_f (∂^2 φ)/〖∂η〗^2
∂φ/∂t=1/S_c (∂^2 φ)/〖∂η〗^2 +S_r (∂^2 θ)/〖∂η〗^2
when η→0
when η→0:u=sin(wt),θ=1 ,φ=1
when η→∞:u→0,θ→0,φ→0
But,i’m confused how to set boundary conditions in this code. Also i want to set these skin friction and nusselt number expressions: C_f=-(∂u/∂η) at η=0
Nu〖Re〗_x^(-1)=-(∂θ/∂η) at η=0
in this code;
function yy= crank
n=105;
h=0.01;
m=0.1;
B=0; %B=0.0001;
Gr=2; %Gr=0.1=Delt. Gs=Dels
A=2; %A=alpha
Pr=5;
Gs=1; %3.5
Sc=1.5;
M=9; %M=30,12
Df=5;
Sr=10;
Ec=0.1; % 0.01-0.1(0-1)
t(1)=0;
s(1)=0;
for j=2:n
t(j)=t(j-1)+h;
end
for j=2:n
s(j)=s(j-1)+m;
end
k=1;
for j=1:n
aa(j,k)=exp(-0.15*Pr^2*t(j)*Ec^-0.1/Gr^0.1); % Guess M*Gr*Gs*S
bb(j,k)=exp(-t(j)*Sc/(Ec)^0.01*M); %bb(j,k)=exp(-Pr^2*t(j)/Ec);
cc(j,k)=exp(-Sc*2*t(j)/Ec^0.3);
end
a{1,k}=[1 0 0;0 1 0;0 0 1];
for j=2:n-1
a{j,k}=[-(1/h)-(1/h^2)-2*A-M Gr Gs;0 -(1/h)-(1/Pr*(h^2)) -Df/h^2;0 -Sr/h^2 -(1/h)-(1/Sc*h^2)];
end
a{n,k}=[1 0 0;0 1 0;0 0 1];
for j=2:n-1
b{j,k}=[(1/2*h^2)+(A/h^3) 0 0;0 1/2*Pr*h^2 Df/2*h^2;0 Sr/2*h^2 1/2*Sc*h^2]; % Lower diagonal entries
end
b{n,k}=[0 0 0;0 0 0;0 0 0];
c{1,k}=[0 0 0;0 0 0;0 0 0];
for j=2:n-1
c{j,k}=[(1/2*h^2)+(A/h^3) 0 0;0 1/2*Pr*h^2 Df/2*h^2;0 Sr/2*h^2 1/2*Sc*h^2];
end
r1(1,k)=0;
r2(1,k)=0;
r3(1,k)=0;
for j=2:n-1
r1(j,k)=-(exp(t(j))/h)-(1/2*h^2)*(1)*(exp(-t(j+1))-2*exp(-t(j))+exp(-t(j-1)))-Gr*exp(-t(j))-Gs*exp(t(j))+M*exp(-t(j))-(A/2*h^3)*(-exp(-t(j+1))+2*exp(-t(j))-exp(-t(j-1)))…
+(aa(j,k)/h)-(1/2*h^2)*(1)*(aa(j+1,k)-2*aa(j,k)+aa(j-1,k))-(A/2*h^3)*( aa(j+1,k)^2-2*aa(j+1,k)*aa(j,k))-Gr*bb(j,k)-Gs*cc(j,k)+M*aa(j,k);
r2(j,k)=(exp((-t(j)))/h)-(1/2*Pr*h^2)*(exp(-t(j+1))-2* exp(-t(j))+exp(-t(j-1)))-Df*(1/2*h^2)*(exp(-t(j+1))-2* exp(-t(j))…
+exp(-t(j-1)))+(bb(j,k)/h)-(1/2*Pr*h^2)*(bb(j+1,k)-2*bb(j,k)+bb(j-1,k))-(Df/2*h^2)*(cc(j+1,k)-2*cc(j,k)+cc(j-1,k));
r3(j,k)=-(exp((-t(j)))/h)-(1/2*Sc*h^2)*(exp(-t(j+1))-2*exp(-t(j))+exp(-t(j-1)))-(Sr/2*h^2)*(exp(-t(j+1))-2*exp(-t(j))+exp(-t(j-1)))+ (cc(j,k)/h)…
-(1/2*Sc*h^2)*(cc(j+1,k)-2*cc(j,k)+cc(j-1,k))-(Sr/2*h^2)*(bb(j+1,k)-2*bb(j,k)+bb(j-1,k));
end
r1(n,k)= 0;
r2(n,k)=0;
r3(n,k)=0;
for j=1:n
rr{j,k}=[r1(j,k);r2(j,k);r3(j,k)];
end
gamma{1,k}=inv(a{1,k})*c{1,k};
for j=2:n-1
a{j,k}=a{j,k}-(b{j,k}*gamma{j-1,k});
gamma{j,k}=inv(a{j,k})*c{j,k};
end
y{1}=inv(a{1})*rr{1};
for j=2:n
y{j}=inv(a{j})*(rr{j}-b{j}*y{j-1});
end
x{n}=y{n};
for j=n-1:-1:1
x{j}=y{j}-(gamma{j})*x{j+1};
end
for j=n:-1:1
u(j,k)=x{j}(1,1);
z(j,k)=x{j}(2,1);
q(j,k)=x{j}(3,1);
end
for j=1:n
xx(j,k)= aa(j,k)+u(j,k);
yy(j,k)= bb(j,k)+z(j,k);
zz(j,k)= cc(j,k)+q(j,k);
end
for j=1:n-1
%U(j)=-((1+(1/B))*((xx(j+1)-xx(j))/h))
%U(j)=-((yy(j+1)-yy(j))/2*h)
end
%ee=meshgrid(xx);
%V=trapz(ee);
%eee=meshgrid(V);
%[XX,YY]=meshgrid(t,s);
%XXX=meshgrid(xx);
%surf(XX,YY,XXX)
%contour(eee)
figure(1)
plot(t,xx,’LineWidth’,2,’MarkerSize’,16,’linestyle’,’-‘,’color’,’k’)
xlabel(‘eta’,’Interpreter’,’tex’,’FontSize’,16,’FontWeight’,’bold’);
ylabel(‘u’,’Interpreter’,’tex’,’FontSize’,16,’FontWeight’,’bold’);
grid on
hold on
figure(2)
plot(t,yy,’LineWidth’,2,’MarkerSize’,16,’linestyle’,’-‘,’color’,’k’)
xlabel(‘eta’,’Interpreter’,’tex’,’FontSize’,16,’FontWeight’,’bold’);
ylabel(‘Theta(eta)’,’Interpreter’,’tex’,’FontSize’,16,’FontWeight’,’bold’);
grid on
hold on
figure(3)
plot(t,zz,’LineWidth’,2,’MarkerSize’,16,’linestyle’,’–‘,’color’,’r’)
xlabel(‘eta’,’Interpreter’,’tex’,’FontSize’,16,’FontWeight’,’bold’);
ylabel(‘Phi(eta)’,’Interpreter’,’tex’,’FontSize’,16,’FontWeight’,’bold’);
grid on
hold on system of pde’s, crank nicloson method in matlab, skin friction and nusselt number MATLAB Answers — New Questions
