how to give the convective boundary condition and guessing function in bvp4c

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);
Error using bvp4c
Unable to solve the collocation equations -- a singular Jacobian encountered.
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

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.
Respected Sir,
Thanks for your response. I have attached the original problem. Kindly make a note on it.
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
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
theta'(0)=-(Bi/A5)*(1-theta(0)), theta(infty)=0
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.

Sign in to comment.

Answers (0)

Categories

Asked:

on 23 Mar 2023

Edited:

on 25 Mar 2023

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!