Why am I getting error using assignin and syms while solving pde using pdepe
Show older comments
So I am trying to solve the non dimensionalized form of an equation which describes pore diffusion resistance combined with surface kinetics. So the equation is a partial differential equation and its boundary condition is a system of coupled ODEs. I have attached the file regarding the equation above. Refer to the non dimensional section for the equations. Here is the code that I have written which uses pdepe and dsolve.
function coupled_pde_ode()
% Define parameters
epsilon_p = input('Porosity (epsilon_p): ');
De = input('Effective diffusivity (De, in m^2/s): ');
R = input('Radius of the pellet (R, in m): ');
kf = input('Mass transfer coefficient (kf, in m/s): ');
kR = input('Reaction rate constant (kR, in 1/s): ');
C0 = input('Initial concentration in pellet pores (C0, arbitrary units): ');
mp = input('Weight of catalyst beads (mp, in kg): ');
rho_p = input('Overall density of catalyst pellets (rho_p, in kg/m^3): ');
V = input('Volume of reactor (V, in m^3): ');
% Derived non-dimensional parameters
phi = (R / 3) * sqrt(kR / De); % Thiele modulus
Bi = kf * R / De; % Biot number
alpha = (3 * epsilon_p * mp) / (V * rho_p); % Coupling parameter
tau_end = 10; % End time in dimensionless units
% Non-dimensionalized spatial and temporal discretization
rbar = linspace(0, 1, 50); % Non-dimensional radius (0 <= rbar <= 1)
tau = linspace(0, tau_end, 100); % Non-dimensional time
% Solve PDE for Cpbar using pdepe
m = 2; % Spherical symmetry
sol_pde = pdepe(m, ...
@(rbar, tau, Cpbar, dCpbar_drbar) pde_eqn(rbar, tau, Cpbar, dCpbar_drbar, epsilon_p, phi), ...
@initial_condition, ...
@(xl, CpbarL, xr, CpbarR, tau) boundary_conditions(xl, CpbarL, xr, CpbarR, tau, Bi, phi), ...
rbar, tau);
% Extract Cpbar at rbar = 1
Cpbar_surface = sol_pde(:, end); % Concentration at rbar = 1 over time
% Solve coupled ODE for Cbbar using dsolve
syms Cbbar(tau)
dCbbar_dtau = -alpha * (Bi * (Cbbar - Cpbar_surface') - 3 * phi^2 * Cpbar_surface');
Cbbar_sol = dsolve(dCbbar_dtau, Cbbar(0) == 1);
% Evaluate Cbbar over time (not plotted)
tau_values = tau; % Time grid
Cbbar_values = double(subs(Cbbar_sol, tau, tau_values));
% Plot results: Concentration vs. Radius
figure;
for i = 1:10:length(tau)
plot(rbar, sol_pde(i, :), 'DisplayName', sprintf('\\tau = %.1f', tau(i)));
hold on;
end
xlabel('\bar{r} (Dimensionless Radius)');
ylabel('\bar{C}_p (Dimensionless Pore Concentration)');
legend;
title('Pore Concentration Profiles (\bar{C}_p) vs. Radius');
grid on;
hold off;
% Nested functions
function [c, f, s] = pde_eqn(rbar, tau, Cpbar, dCpbar_drbar, epsilon_p, phi)
% Non-dimensional PDE coefficients
c = epsilon_p; % Time-dependent term coefficient
f = dCpbar_drbar; % Diffusion term coefficient
s = -9 * phi^2 * Cpbar; % Source term (reaction rate)
end
function Cpbar0 = initial_condition(rbar)
% Non-dimensional initial condition
Cpbar0 = ones(size(rbar)); % Initial concentration equals C0 (normalized to 1)
end
function [pl, ql, pr, qr] = boundary_conditions(xl, CpbarL, xr, CpbarR, tau, Bi, phi)
% Boundary conditions for the non-dimensional problem
% At rbar = 0 (center of the sphere), Neumann BC: dCpbar/drbar = 0
pl = 0;
ql = 1;
% At rbar = 1 (outer surface of the sphere), flux term
pr = Bi * (1 - CpbarR) - 3 * phi^2 * CpbarR;
qr = 1;
end
end
Here are the values I entered
Porosity (epsilon_p):
0.5
Effective diffusivity (De, in m^2/s):
1.5e-10
Radius of the pellet (R, in m):
1e-3
Mass transfer coefficient (kf, in m/s):
2.8e-3
Reaction rate constant (kR, in 1/s):
5e-12
Initial concentration in pellet pores (C0, arbitrary units):
1
Weight of catalyst beads (mp, in kg):
0.1e-3
Overall density of catalyst pellets (rho_p, in kg/m^3):
850
Volume of reactor (V, in m^3):
0.1
Here is the error message I'm receiving
Error using assignin
Attempt to add "Cbbar" to a static workspace.
Error in syms (line 331)
assignin('caller', name, xsym);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in Testing (line 35)
syms Cbbar(tau)
Accepted Answer
Here is the code for the dimensional model.
It should be no problem to rewrite it in the non-dimensional form.
Note that the sign in S1c of your attached document is incorrect:
it should be
kf*(Cp-Cb)
instead of
kf*(Cb-Cp)
% Define parameters
epsilon_p = 0.5;%input('Porosity (epsilon_p): ');
De = 1.5e-10;%input('Effective diffusivity (De, in m^2/s): ');
R = 1e-3;%input('Radius of the pellet (R, in m): ');
kf = 2.8e-3;%input('Mass transfer coefficient (kf, in m/s): ');
kR = 5.0e-12;%input('Reaction rate constant (kR, in 1/s): ');
C0 = 1;%input('Initial concentration in pellet pores (C0, arbitrary units): ');
mp = 1e-4;%input('Weight of catalyst beads (mp, in kg): ');
rho_p = 850;%input('Overall density of catalyst pellets (rho_p, in kg/m^3): ');
V = 1e-1;%input('Volume of reactor (V, in m^3): ');
% spatial and temporal discretization
nr = 100;
r = linspace(0, R, nr).'; % radius (0 <= rbar <= 1)
dr = r(2) - r(1);
tend = 100; % End time
nt = 10;
tspan = linspace(0, tend, nt); % time
% Initial conditions
Cp0 = zeros(nr,1);
Cb0 = 1.0;
y0 = [Cp0;Cb0];
[T,Y] = ode15s(@(t,y)fun(t,y,nr,r,dr,epsilon_p,De,R,kf,kR,C0,mp,rho_p,V),tspan,y0);
Cp = Y(:,1:end-1);
Cb = Y(:,end);
figure(1)
plot(r,Cp)
figure(2)
plot(T,Cb)
function dydt = fun(t,y,nr,r,dr,epsilon_p,De,R,kf,kR,C0,mp,rho_p,V)
Cp = y(1:end-1);
Cb = y(end);
dCp_dt = zeros(nr,1);
dCp_dt(1) = 3/((r(1)+dr/2)^3-r(1)^3)*De*...
((r(1)+dr/2)^2*(Cp(2)-Cp(1))/dr)-...
kR * Cp(1);
dCp_dt(2:nr-1) = 3./((r(2:nr-1)+dr/2).^3-(r(2:nr-1)-dr/2).^3)*De.*...
((r(2:nr-1)+dr/2).^2.*(Cp(3:nr)-Cp(2:nr-1))/dr-...
(r(2:nr-1)-dr/2).^2.*(Cp(2:nr-1)-Cp(1:nr-2))/dr)-...
kR * Cp(2:nr-1);
dCp_drR = -(kf*(Cp(nr)-Cb) - R/3*kR*Cp(nr))/De;
dCp_dt(nr) = 3/(r(nr)^3-(r(nr)-dr/2)^3)*De*...
(r(nr)^2*dCp_drR - ...
(r(nr)-dr/2)^2*(Cp(nr)-Cp(nr-1))/dr)-...
kR * Cp(nr);
dCp_dt = dCp_dt/epsilon_p;
dCb_dt = -1/V * mp/rho_p * 3/R * De * dCp_drR;
dydt = [dCp_dt;dCb_dt];
end
6 Comments
Aditya
on 16 Dec 2024
In order to solve (S4) over time, you should insert (S3c) into (S4). (As mentionned above, there is a sign error in (S3c) in the convective flux term.)
Then, from the solution of "pdepe", you have to extract the values for Cpbar at rbar = 1 over time and use interpolation in ode_rhs to the times "tau" requested by ode45 (you shouldn't name them "tau" in the list of input parameters as they differ from the tau-vector used as output times for "pdepe").
So all in all: you might succeed, but iteration usually is the last refuge. The cleanest way is to define the PDE and ODE together and use an ODE integrator that solves all the resulting equations simultaneously.
I think that
dCpbar_drbar_guess = (Cpbar_new(end) - Cpbar_new(end-1)) / (rbar(end) - rbar(end-1));
is incorrect since
Cpbar_new(end) - Cpbar_new(end-1)
looks like a difference in time, divided by
rbar(end) - rbar(end-1)
which is a difference in space.
Walter Roberson
on 16 Dec 2024
You need
xlabel('\bar{r} (Dimensionless Radius)', 'interpreter', 'latex');
ylabel('\bar{C}_p (Dimensionless Pore Concentration)', 'interpreter', 'latex');
title('Final Pore Concentration Profile (\bar{C}_p) vs Radius', 'interpreter', 'latex');
If kR*Cp in your document must be replaced by kR*Cp.^n for an n-th order reaction, the changes you made should be correct.
For the discretization, see Chapter 3 of the enclosed document.
Aditya
on 18 Dec 2024
More Answers (1)
tau is a vector of values in your code:
tau = linspace(0, tau_end, 100); % Non-dimensional time
You can't use it as a symbolic variable afterwards without redefining it as symbolic:
syms tau Cbbar(tau)
But note that you now have overwritten the original tau.
Further, replace
dCbbar_dtau = -alpha * (Bi * (Cbbar - Cpbar_surface') - 3 * phi^2 * Cpbar_surface');
by
dCbbar_dtau = diff(Cbbar,tau) == -alpha * (Bi * (Cbbar - Cpbar_surface') - 3 * phi^2 * Cpbar_surface');
to define your differential equation.
But this will only work for Cpbar_surface being a single value, not a vector of values that depend on time.
Thus you won't succeed to integrate
dCbbar_dtau = -alpha * (Bi * (Cbbar - Cpbar_surface') - 3 * phi^2 * Cpbar_surface');
with a time-dependent function Cpbar_surface(t) (in this case a vector of values) symbolically.
Use pde1dm if you want to solve a PDE with a coupled ODE:
Or discretize your PDE in space, add the ODE and solve the complete system of ordinary differential equations using "ode15s". Look up "method-of-lines" for more details.
Categories
Find more on Numerical Integration and Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
