The convergence problem you are seeing is due to the strong boundary layer effect
near 
, i.e. the solution is undergoing a sharp change in this region. The
, i.e. the solution is undergoing a sharp change in this region. Thediscretization method in pdepe requires a fine mesh in this region to obtain a
satisfactory solution.
If you solve this system for a large value of Ra, and examine the solution closely in this
region, you will see some non-physical oscillations in the solution values. If you increase
Ra still further, these spurious oscillations become larger, eventually leading to the
convergence failure you are seeing.
Although a very fine mesh is needed near , a much coarser mesh is acceptable
, a much coarser mesh is acceptableelsewhere. In my version of your example, which I show below, I have exponentially
graded the mesh to produce a much denser mesh around .
.function matlabAnswers_10_11_2019
%% In this code u1= v (vertical velocity), u2= Temperature, y is the vertical dimension
%% setting up the problem parameters
L=2; %% total depth, in positive units
n=2000;
if(0)
y= linspace (0,L,n); %% space mesh
else
y=expMesh(L, n);
end
t= linspace(0,0.00015,3); %% time mesh
m=0; %% part of pdepe arguments, we have cartisian coordinate system, otherwise m changes
tic
sol = pdepe(m,@heatpde,@heatic,@heatbc,y,t); %% solving system of PDEs
toc
%% for mapping the result
u1 = sol(:,:,1); %% obtained Velocity profile
u2 = sol(:,:,2); %% obtained Temperature Profile
%figure; plot(u2,-y); grid;
figure; plot(u2(end,:),-y); grid;
%plot(u1,-y)
end
%% defining PDE function
function [c,f,s] = heatpde(y,t,u,dudy)
%properties of fluid under consideration
Ra= 1e7; %% Raynold Number
Ra = 1e5;
Pr= 7; %% Prandlt Number
R= 1; %% density
G= 10; % gravitational acceleration
% function matrices
c=[1;1];
f= [Pr;1].*dudy; %% flux term
F1= Ra *Pr *u(2);
F2= Ra *Pr * R * G;
F3= u(1)* dudy(1);
F4= 2.718^(-y);
F5= u(1)* dudy(2);
s= [(F1+F2-F3);(F4-F5)]; %% source term
end
%% Initial Conditions
function u0= heatic(y)
u0 = [0;10]; %%initial velocity and temperature
end
%% Boundary Conditions
function [pt,qt,pb,qb] = heatbc (yt,ut,yb,ub,t)
F5= -2.718^(-2);
pt= [0; 0]; %%top B.C,
qt= [1; 1]; %% top B.C,
pb= [ub(1); F5]; %%bottom B.C,
qb= [0; 1]; %% bottom B.C, p+(q*f)=0
end
function x=expMesh(L, n)
x=linspace(0,L, n);
h=exp(-x);
hp=L/sum(h)*h;
x=[0 cumsum(hp)];
%figure; plot(1:n+1, x, 'o');
%z = zeros(n+1,1);
%figure; plot(x, z, 'o');
end
