6 views (last 30 days)

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

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

Opportunities for recent engineering grads.

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

Start Hunting!
## 0 Comments

Sign in to comment.