# FEM for a biphasic drying problem

I have the following PDE originating from a biphasic drying model, where $$\xi \in [0,\Xi]$$ is a radial coordinate attached to the dry skeleton of a wet cylindrical body: $$\frac{\partial u}{\partial t} = \frac{1}{\xi}\frac{\partial}{\partial\xi}\left(\xi D_\xi \frac{\partial u}{\partial \xi}\right), \quad D_\xi = \frac{D}{(1+\alpha u)^2} \left(\frac{r}{\xi}\right)^2$$ The scalar $$u = u(\xi,t) = \rho_\textrm{l}(r,t)/\rho_\textrm{s}(r,t)$$ is the moisture content on a dry basis, defined as the ratio of volumetric concentration of liquid and solid, $$D$$ is an effective diffusivity, $$\alpha$$ is a shrinkage factor, and $$r = r(u(\xi),\xi)$$, $$r \in [0,R(t)]$$, is the radial position corresponding to $$\xi$$ in the wet reference frame (i.e. the Eulerian or geometric reference frame). At $$t = 0$$ we assume a homogeneous moisture distribution $$u(\xi,0) = u_0$$.

At $$\xi = r = 0$$, the moisture distribution is symmetric: $$\frac{\partial u}{\partial \xi} = 0.$$ At the surface $$\xi = \Xi$$ (or $$r = R(t)$$) either a Dirichlet or Robin condition is prescribed depending on the magnitude of the Biot number. Both are captured by the equation: $$-\rho_\textrm{s}^* D_\xi \frac{\partial u}{\partial \xi} = k_Y(Y^*(u) - Y_\textrm{air}).$$ When ventilation in the dryer is strong, the Biot number, $$\textit{Bi} = k_Y R(t)/D$$ will approach infinity, leaving the equilibrium relation $$Y^*(u) = Y_\textrm{air}$$ which can be solved for $$u|_\Xi = u^*$$. The function $$Y^*$$ is a constitutive equation also known as the moisture sorption isotherm.

The relation between $$r$$ and $$\xi$$ is found from the mass balance of solid over a differential element (an annulus), $$dm_\textrm{s} = \textrm{constant} \quad\longrightarrow\quad \rho_\textrm{s} (2\pi r) dr = \rho_{\textrm{s}}^* (2\pi \xi)d\xi,$$ where $$\rho_\textrm{s}$$ is the concentration of solids, and $$\rho_\textrm{s}^*$$ is the partial density of pure solids. In an ideally-shrinking biphasic system, the two quantities are related by a constitutive equation: $$\frac{\rho_\textrm{s}}{\rho_\textrm{s}^*} = \frac{1}{1 + \alpha u}$$

By taking the integral over the differential element, we can find the wet radius: $$r^2 = \xi^2 + 2\alpha\int_0^\xi u(\xi') \xi' d\xi' \longrightarrow \left(\frac{r}{\xi}\right)^2 = 1 + \alpha \frac{2}{\xi^2} \int_0^\xi u(\xi') \xi' d\xi'$$

To solve the model I'd like to use the method of lines in combination with a finite-element discretization of the governing equation. It appears like at each time step I need to compute the cumulative integral of the moisture content to obtain the values $$r(\xi)$$. (For example the MATLAB cumtrapz function could be used for this purpose.)

Alternatively, taking the derivative of $$r$$ with respect to time, and the fact that $$\partial \xi /\partial t = 0$$, we can do the following manipulation: $$2r\frac{\partial r}{\partial t} = 2\alpha \int_0^\xi \frac{\partial u(\xi',t)}{\partial t}\xi'd\xi' = 2\alpha \xi D_\xi \frac{\partial u}{\partial \xi},$$ giving, $$\frac{\partial r}{\partial t} = \alpha \frac{\xi}{r} D_\xi \frac{\partial u}{\partial \xi}.$$ This equation gives the trajectories of $$r(\xi,t)$$ and lends itself to discretization using the method of lines as an additional set of ODEs for the position of the finite element edges. With a $$C^0$$ FEM basis however, the derivative $$\partial_\xi u$$ is discontinuous between elements, hence it would appear this method is not applicable.

What options do I have to circumvent this problem and which of them would you recommend?

My current list of options is:

• Use the integral relation and forget about $$\partial r/\partial t$$
• Use an $$L_2$$-projection of the derivative $$\partial_\xi u$$ onto the original $$C^0$$ basis
• Use a $$C^1$$ basis, where the derivatives at the edges are continuous
• Use FDM or FVM instead

Any other thoughts will be also appreciated.

For better imagination of what the solutions tend to look like (here shown for a constant Dirichlet boundary condition $$u^* = 0.1$$, and initial concentration $$u_0 = 0.3$$):

Note the movement of the surface in the bottom right corner.

• In your first equation, what is x? Oct 25, 2023 at 13:16
• It was a typo. Sorry about that. Oct 25, 2023 at 13:28
• I don't see any reason why either of your two proposed approaches won't work. It is not uncommon in FEM for spatial derivatives of a dependent variable to be discontinuous at the element nodes; this is one of the strengths of FEM. You typically need to compute these derivatives only at interior points in an element. Oct 26, 2023 at 13:36
• Can you include the boundary and initial conditions for your problem? Oct 26, 2023 at 13:36
• I've added the BC's. Are you implying I only need to evaluate $r$ at the quadrature points, so the edge values are not needed? While I tend to solve for $u$ as a function of $\xi$, for plotting it's convenient to know the $r$ values at the cell edges. Oct 26, 2023 at 14:45

I solved your equations using the pdepe function in MATLAB and your second approach. That is, I included the second PDE that describes the evolution of $$r$$ with time. This two-PDE approach is more straightforward to implement with pdepe. The pdepe function basically uses a low-order finite element approximation. As I said previously, I believe the choice between your two approaches is mostly due to computational convenience.

This analysis produces the figure below that is comparable to yours.

Below is the MATLAB script for this problem.

function cse_10_26_2023_pdepe
m = 1;
nx=100;
r0=1e-3;
Xi=0.00083;
xi = linspace(0,Xi,nx);
t=[0 30 60 120 240 480 960 1600 3200 6400];
nt=length(t);
D=1e-10;
alpha=1.5;
u0=.3;

pdef=@(x,t,u,DuDx) pdefun(x,t,u,DuDx,D,alpha);
bcf=@(xl,ul,xr,ur,t) bcfun(xl,ul,xr,ur,t);
icf=@(xi) icfun(xi,u0,Xi,r0);

opts=odeset('abstol', 1e-12);
sol = pdepe(m,pdef,icf,bcf,xi,t,opts);
u=sol(:,:,1);
r=sol(:,:,2);

figure; hold on;
for i=1:nt
plot(r(i,:), u(i,:));
end
hold off;
title 'u at different times';
axis([0 r0 .1 .35]);
timeLabels=arrayfun(@num2str, round(t), 'UniformOutput', 0);
legend(timeLabels);
xlabel 'r';

figure; plot(xi, r(end,:)); title 'r at final time';
figure; plot(t, u(:,1)); title 'u at center as a function of time';
figure; plot(t, r(:,end)); title 'r at outer surface as a function of time';
disp(r(:,end)');
end

function [c,f,s] = pdefun(xi,t,U,DUDxi,D,alpha)
u=U(1); r=U(2);
DuDxi=DUDxi(1);
Dxi=dxi(D, alpha, xi, r, u);
c = [1 1]';
f = [Dxi*DuDxi 0]';
s = [0 alpha*xi/r*Dxi*DuDxi]';
end

function U0 = icfun(xi,u0,Xi,r0)
U0=[u0 xi/Xi*r0]';
end

function [pl,ql,pr,qr] = bcfun(rCenter,uCenter,rOuter,uOuter,t)
pl = [0 uCenter(2)]';
ql = [1 0]';
pr = [uOuter(1)-.1 0]';
qr = [0 1]';
end

function Dxi=dxi(D, alpha, xi, r, u)
Dxi=D/(1+alpha*u)^2*(r/xi)^2;
end

• Thanks for the wonderful demonstration. The pdepe routine uses a Petrov-Galerkin method when the domain includes the singular point $r = 0$. Your boundary functions for the right side sets $f(\Xi,t,u,\partial_\xi u) = 0$ (zero-flux), however f = 0 from pdefun. This strikes me as slightly odd, even though "do nothing" is precisely what I'd want for $r$ at that boundary. Is setting qr = 1; pr = 0 always granted to give behavior consistent with a natural boundary condition? Oct 29, 2023 at 17:55
• Yes, setting the BC that way is simply a trick to satisfy pdepe requirements. Oct 29, 2023 at 18:15
• My full problem also has a coupled ODE (a lumped energy equation for the temperature), but it appears like pdepe can't be used in this case. Would the pde1dm solver work in this case? Feb 5 at 16:22
• I'm not sure what you mean by a lumped energy equation but pde1dm does allow for coupled ODE. A common use for this is to account for a non-standard BC. If you run into problems, post a mathematical description of your equations along with your matlab code. Feb 5 at 21:29