how to give the convective boundary condition and guessing function in bvp4c
조회 수: 3 (최근 30일)
이전 댓글 표시
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
답변 (0개)
참고 항목
카테고리
Help Center 및 File 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!