MATLAB Answers

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

6 views (last 30 days)
Muhammad Taha Manzoor
Muhammad Taha Manzoor on 10 Oct 2019
Edited: Bill Greene on 12 Oct 2019
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

  0 Comments

Sign in to comment.

Answers (1)

Bill Greene
Bill Greene on 12 Oct 2019
Edited: Bill Greene on 12 Oct 2019
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

  0 Comments

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!