how to give the convective boundary condition and guessing function in bvp4c
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 23 Mar 2023
We must see a mathematical description of your problem (equations, boundary conditions).
If we tried to extract this information from your code, we also extracted the possible errors you made.
N Nithya
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)
Categories
Find more on Polar Plots 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!