PDE Toolbox: Solver produces error depending on the parameters and initial conditions

Question

Anselm Dreher 2021-10-5

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1466991-pde-toolbox-solver-produces-error-depending-on-the-parameters-and-initial-conditions

编辑： Anselm Dreher 2021-10-19

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

C0 = Cfeed * 2e-1

and above, but not for

C0 = Cfeed * 1e-1

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

% Create PDE model with 1 equation

model = createpde(1);

L = 1; % length

H = 1e-1; % arbitrary height

D = 1e-2; % Diffusion coefficient m^2/s

v = 4e-6; % velocity m/s

nfeed = 0.1; % Feed mole flow mol/s

n0 = nfeed * 1e-1; % Initial condition mol/s;

kinParam = 5e-2; % kinetic parameter

hmax = H; % Size of FEM element

% Calculate Peclet-Number, ratio of convection and diffusion

Pe = v * hmax / (2*D);

% Time vector

tspan = linspace(0,100,100);

% define geometry for fluid flow

F_Rct = [3 4 0 L L 0 0 0 H H]'; %Description Matrix

gd = F_Rct;

sf = 'F1';

ns = ('F1')';

geo = decsg(gd,sf,ns);

% assign geometry to model

geometryFromEdges(model,geo);

% plot geometry

figure(1)

pdegplot(model,'EdgeLabels','on','FaceLabels','on')

% generate mesh

generateMesh(model,'Hmax',hmax);

pdemesh(model)

% assign coefficients of the PDE

% for coefficient d assign 1 for non stationary problem

% for coefficient f pass the additional parameter v for fluid velocity

specifyCoefficients(model,'m',0,...

'd',1,...

'c',D,...

'a',0,...

'f',@(location,state)f_coeff (location,state,v,kinParam));

% define boundary conditions, Edge 4 is inlet

applyBoundaryCondition(model,'dirichlet','Edge',4,'u',nfeed);

applyBoundaryCondition(model,'neumann','Edge',[1 2 3]);

% define initial conditions

setInitialConditions(model,n0);

% solve PDE

results = solvepde(model,tspan);

% Plot solution

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 defining coefficient f

function f = f_coeff(~,state,v,kinParam)

% convective term state.ux devided by arbitrary factor state.u

% followed by source term (reaction)

f = -v * state.ux ./ state.u - kinParam * state.u;

end

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Anselm Dreher 2021-10-19

编辑：Anselm Dreher 2021-10-19

One possible workaround is to modify it to a stationary problem which makes the script working properly. just set parameter d to 0 and remove tspan from solvepde. with different parameters (that dont work with the non-stationary code) one even gets a meaningful concentration profile

Parameters:

L = 1; % length

H = 1e-1; % arbitrary height

D = 5e-8; % Diffusion coefficient m^2/s

v = 3e-6; % velocity m/s

Cfeed = 0.1; % Feed concentration mol/m^3

C0 = Cfeed * 1e-1; % Initial condition mol/m^3;

kinParam = 5e-4;

plot:

figure(2)

pdeplot(model,'XYData',results.NodalSolution(:,1))

caxis([0 0.1])