PDEPE: Unable to meet integration tolerances for large values of Ra

Hi,

The following code runs fine when my Ra values are low or I reduce the t to very small values. But as soon as I incerese t or Ra the code gives the following error.

"Warning: Failure at t=5.652476e-05. Unable to meet integration tolerances without reducing the step size below the smallest

value allowed (1.084202e-19) at time t. "

I am trying to solve convective heat transfer equation in a vessel heated from the bottom. B.C are that at the top I have a stress free and adiabatic surface. At the bottom I have a constant flux and no-slip.

%% 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

y= linspace (0,L,30); %% space mesh

t= linspace(0,0.00015,3); %% time mesh

m=0; %% part of pdepe arguments, we have cartisian coordinate system, otherwise m changes

sol = pdepe(m,@heatpde,@heatic,@heatbc,y,t); %% solving system of PDEs

%% for mapping the result

u1 = sol(:,:,1); %% obtained Velocity profile

u2 = sol(:,:,2); %% obtained Temperature Profile

plot(u2,-y)

%plot(u1,-y)

%% defining PDE function

function [c,f,s] = heatpde(y,t,u,dudy)

%properties of fluid under consideration

Ra= 1*10^7; %% Raynold Number

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

Answer 1

[Bill Greene]

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

discretization 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

elsewhere. 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