% this function estimates the correct value of a for big Pr cases
function f=big_Pr(a)
x0=[0 0 a(1) 1 a(2)];
[t,Y]=ode15s(@governing_equation,[0 7],x0);
% check the BC's at infinity approach zero
f(1)=Y(end,2);
f(2)=Y(end,4)^2+Y(end,5)^2;
end
% this function defines the similarity equations as defined by
% Bejan using correct scaling and Ostrach using incorrect scaling
function xdot=governing_equation(t,x)
global Pr
% Based on Bejan's scaling similarity equation
% The state-space derivation is shown in the document
xdot(1)=x(2);
xdot(2)=x(3);
xdot(3)=-(1/Pr)*(0.5*x(2)^2-(3*x(1)*x(3))/4)+x(4);
xdot(4)=x(5);
xdot(5)=3*x(1)*x(5)/4;
xdot=xdot';
% Based on Ostrach's scaling similarity equation
% % xdot(1)=x(2);
% % xdot(2)=x(3);
% % xdot(3)=-3*x(1)*x(3)+2*x(2)^2-x(4);
% % xdot(4)=x(5);
% % xdot(5)=-3*Pr*x(1)*x(5);
% % xdot=xdot';
end
global Pr
Pr = 1
options=optimset('TolFun',1e-6,'TolX',1e-6);
a=fsolve(@big_Pr,[0.5 -0.6],options);
a
x0=[0 0 a(1) 1 a(2)];
[t,Y]=ode15s(@governing_equation,[0 4],x0);
figure(1);
plot(t,-Y(:,2),'r')
hold on
figure(2);
plot(t,Y(:,4),'r')
hold on
Pr=100
options=optimset('TolFun',1e-6,'TolX',1e-6);
a=fsolve(@big_Pr,[0.25 -0.7],options);
a
x0=[0 0 a(1) 1 a(2)];
[t,Y]=ode15s(@governing_equation,[0 4],x0);
figure(1);
plot(t,-Y(:,2),'b')
figure(2);
plot(t,Y(:,4),'b')
Pr=1000
options=optimset('TolFun',1e-6,'TolX',1e-6);
a=fsolve(@big_Pr,[0.1 -0.65],options);
a
x0=[0 0 a(1) 1 a(2)];
[t,Y]=ode15s(@governing_equation,[0 4],x0);
figure(1);
plot(t,-Y(:,2),'k')
figure(2);
plot(t,Y(:,4),'k')
Pr=10
options=optimset('TolFun',1e-6,'TolX',1e-6);
a=fsolve(@big_Pr,[0.5 -0.52],options);
a
x0=[0 0 a(1) 1 a(2)];
[t,Y]=ode15s(@governing_equation,[0 4],x0);
figure(1);
plot(t,-Y(:,2),'g')
figure(2);
plot(t,Y(:,4),'g')
Pr=0.01
options=optimset('TolFun',1e-6,'TolX',1e-6);
a=fsolve(@small_Pr,[0.15 -0.17],options);
a
x0=[0 0 a(1) 1 a(2)];
[t,Y]=ode15s(@governing_equation,[0 4],x0);
figure(1);
plot(t,-Y(:,2),'y')
legend('Pr = 1','Pr = 100','Pr = 1000','Pr = 10','Pr = 0.01')
title('Vertical Velocity Profiles')
xlabel('η = (x/y) Ra_y^(^0^.^2^5^)')
ylabel('G = v(y/α) Ra_y^(^0^.^2^5^)')
figure(2);
plot(t,Y(:,4),'y')
legend('Pr = 1','Pr = 100','Pr = 1000','Pr = 10','Pr = 0.01')
title('Temperature Profiles')
xlabel('η = (x/y) Ra_y^(^0^.^2^5^)')
ylabel('ϴ')
function f=small_Pr(a)
x0=[0 0 a(1) 1 a(2)];
[t,Y]=ode15s(@governing_equation,[0 30],x0);
f(1)=Y(end,2);
f(2)=Y(end,4)^2+Y(end,5)^2;
end
%% I want to run this code, please assemble it so that it will run in 2023b version