how to give the convective boundary condition and guessing function in bvp4c
3 views (last 30 days)
Show older comments
Respected sir,
The following are the problem of mixed convection flow of fluid. I have a mention my code and figure using bvp4c. I am getting velocity profile, but not getting proper temperature profile.kindly give me the suggestion for how to choose the guessing function.
global phi Bi fw pr A M lamda Ri Rd Ec Hg ks kf rowf rows betas betaf sigmaf sigmas rownf cpf cpnf betanf rowcpnf rowcpf rowbetanf rowbetaf A1 A2 A3 A4 A5 A6
%aluminium oxide
phi = 0;
Bi=0.1;
fw=0.1;
pr = 6.2;
A = 0.3;
M = 1.5;
lamda = 0.3;
Ri = 0.5;
Rd = 0.4;
Ec = 0.5;
Hg = -0.4;
ks = 40;
kf = 0.643;
rowf = 997.1;
rows = 3970;
betas = 0.85*10.^-5;
betaf = 21*10.^-5;
sigmaf = 0.05;
sigmas = 1*10.^-10;
rownf = (1-phi)*rowf+phi*rows;
cpnf = 4182;
rowf = 997.1;
cpf = 4179;
betanf = (1-phi)*betaf+phi*betas;
rowcpnf = rownf*cpnf;
rowcpf = rowf*cpf;
rowbetanf = rownf*betanf;
rowbetaf = rowf * betaf;
%knf/kf = 1;
A1 = (1-phi).^2.5;
A2 = rowf/rownf;
A3 = 1+((3*((sigmas/sigmaf)-1)*phi))/(((sigmas/sigmaf)+2)-phi*((sigmas/sigmaf)-1));
A4 = rowbetanf/rowbetaf;
A5 = ((1-phi)+2*phi*(ks/(ks-kf))*log((ks+kf)/2*kf))/((1-phi)+2*phi*(kf/(ks-kf))*log((ks+kf)/2*kf));
A6 = rowcpnf/rowcpf;
%putting the original formula for A5
xmesh = linspace(0,7,20);
solinit1 = bvpinit(xmesh, @guess);
sol1 = bvp4c(@flow,@bvp,solinit1);
y=deval(sol1,xmesh);
figure(1)
plot(sol1.x,sol1.y(4,:),'linewidth',1.5);
%plot(sol.x,sol.y(5,:),'--','linewidth',1.5);
hold on
%constants
function dy = flow(t,y)
global phi Bi fw pr A M lamda Ri Rd Ec Hg ks kf rowf rows betas betaf sigmaf sigmas rownf cpf cpnf betanf rowcpnf rowcpf rowbetanf rowbetaf A1 A2 A3 A4 A5 A6
% for clearer code
F0 = y(1);
F1 = y(2);
F2 = y(3);
G0 = y(4);
G1 = y(5);
%Define a function
%A1*A2*f'''(t)+f(t)*f''(t)-f'(t)^2-A[f'(t)+t/2*f''(t)]-A1*A2*lambda*f'(t)-A2*A3*M*f'(t)+A2*A4*Ri*theta(t)=0
dy(1)= F1;
dy(2)= F2;
dy(3)=( F1.^2-F0*F2+A*[F1+(t/2)*F2]+A1*A2*lamda*F1+A2*A3*M*F2-A2*A4*Ri*G0)/(A1*A2);
%A5*A6*(1/Pr)*theta''(t)+f(t)*theta'(t)-f'(t)*theta(t)-A*[theta(t)+(t/2)*theta'(t)]+A3*A6*M*Ec*f'(t).^2+A6*Hg*theta(t)=0
f(0)=fw; f'(0)=1, f'(infty)=0
dy(4)= G0;
dy(5)=(pr* (F1*G0-F0*G1+A*[G0+(t/2)*G1]-A1*A6*Ec*F2.^2-A3*A6*M*Ec*F1.^2-A6*Hg*G0))/(A6*(A5+(4/3)*Rd));
dy = [dy(1);dy(2);dy(3);dy(4);dy(5)];
end
function res = bvp(ya,yb)
global phi Bi fw pr A M lamda Ri Rd Ec Hg ks kf rowf rows betas betaf sigmaf sigmas rownf cpf cpnf betanf rowcpnf rowcpf rowbetanf rowbetaf A1 A2 A3 A4 A5 A6
res = [ ya(1)- fw
ya(2)-1
% ya(5)+(Bi/A5)*(1-ya(4))
ya(4)-1
yb(2)
yb(4)];
end
function g = guess(t)
global phi Bi fw pr A M lamda Ri Rd Ec Hg ks kf rowf rows betas betaf sigmaf sigmas rownf cpf cpnf betanf rowcpnf rowcpf rowbetanf rowbetaf A1 A2 A3 A4 A5 A6
g=[fw
1
0.0001
exp(-4*t/2)
0.111];
end
3 Comments
Torsten
on 25 Mar 2023
Edited: Torsten
on 25 Mar 2023
I don't see that the equation you post above for theta equals the one you define in your code. The division by (A5*A6*(1/Pr)) is wrong and the term A1*A6*Ec*F2.^2 cannot be found.
Further, dividing the first equation by A1*A2 means /(A1*A2), not /A1*A2.
Answers (0)
See Also
Categories
Find more on Numerical Integration and Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!