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

조회 수: 3 (최근 30일)
N Nithya
N Nithya 2023년 3월 23일
편집: Torsten 2023년 3월 25일
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
N Nithya
N Nithya 2023년 3월 25일
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
Torsten
Torsten 2023년 3월 25일
편집: Torsten 2023년 3월 25일
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.

댓글을 달려면 로그인하십시오.

답변 (0개)

카테고리

Help CenterFile Exchange에서 Numerical Integration and Differential Equations에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by