Hello,
I am trying to solve a PDE system describing a convection - diffusion - problem (i.e. a chemical reactor).
I created a minimum working example to demonstrate my Problem. The PDE describes a 1D reactor with a substance flowing into the reactor with a feed concentration and being transportet by convection and diffusion.
The coefficients of my PDE system require to calculate material properties depending on the current state of the system, for example by dividing by "state.u". Depending on the set parameters for the convective term, the diffusive term and the initial conditions the following warning and error message is produced:
Warning: Failure at t=1.515228e-08. Unable to meet integration tolerances without reducing the step size
below the smallest value allowed (5.293956e-23) at time t.
> In ode15s (line 655)
In pde.EquationModel/solveTimeDependent (line 101)
In pde.PDEModel/solvepde (line 57)
In execute_minEx_error_value_too_small (line 65)
Error using pde.EquationModel/solveTimeDependent (line 103)
Solution failed to reach the requested end time.
Error in pde.PDEModel/solvepde (line 57)
[u,dudt] = self.solveTimeDependent(coefstruct, u0, ut0, tlist, ...
Error in execute_minEx_error_value_too_small (line 65)
results = solvepde(model,tspan);
For the given values of v and D the code works fine for an inital condition of
and above, but not for
and below. It seems to be the case that by reducing the Peclet number (decreasing v, increasing D) the initial conditions can be set to a lower value. Removing the division by state.u in the f coefficient makes the program to be able to run with any C0 value.
Can you help me to understand why the error is produced and how to make my model run for a wider range of parameters? Unfortunately I cannot simulate my system if I cant choose the parameters more freely.
clearvars, close all
clc
model = createpde(1);
L = 1;
H = 1e-1;
D = 1e-2;
v = 4e-6;
nfeed = 0.1;
n0 = nfeed * 1e-1;
kinParam = 5e-2;
hmax = H;
Pe = v * hmax / (2*D);
tspan = linspace(0,100,100);
F_Rct = [3 4 0 L L 0 0 0 H H]';
gd = F_Rct;
sf = 'F1';
ns = ('F1')';
geo = decsg(gd,sf,ns);
geometryFromEdges(model,geo);
figure(1)
pdegplot(model,'EdgeLabels','on','FaceLabels','on')
generateMesh(model,'Hmax',hmax);
pdemesh(model)
specifyCoefficients(model,'m',0,...
'd',1,...
'c',D,...
'a',0,...
'f',@(location,state)f_coeff (location,state,v,kinParam));
applyBoundaryCondition(model,'dirichlet','Edge',4,'u',nfeed);
applyBoundaryCondition(model,'neumann','Edge',[1 2 3]);
setInitialConditions(model,n0);
results = solvepde(model,tspan);
figure(2)
for i = 1:length(tspan)
pdeplot(model,'XYData',results.NodalSolution(:,i))
caxis([0 0.1])
title(['Amount of substance n (' num2str(i) '/' num2str(length(tspan)) ')'])
pause(0.001)
end
function f = f_coeff(~,state,v,kinParam)
f = -v * state.ux ./ state.u - kinParam * state.u;
end