Robin boundary conditions or mixed Dirichlet (prescribed value) and Neumann (flux) conditions are a third type of boundary condition that for example can be used to implement convective heat transfer and electromagnetic impedance boundary conditions. In the following it will be discussed how mixed Robin conditions are implemented and treated in FEATool with an illustrative example (in short Robin conditions can be entered in FEATool as the usual Neumann flux coefficients).
Consider the following one-dimensional convection, diffusion, and reaction model example
\[ \frac{\partial^2 u}{\partial x^2} + u\frac{\partial u}{\partial x} - u = e^{2x} \]
with the two boundary conditions
\[
\left\{\begin{aligned}
u+\frac{\partial u}{\partial x} = 2 &, \qquad x = 0\\\
u = e &, \qquad x = 1
\end{aligned}\right.
\]
and the exact solution $u(x) = e^x$.
The first boundary condition $u+\frac{\partial u}{\partial x} = 2$ at $x=0$ is of the mixed Robin type, prescribing a combination of both the solution value (Dirichlet) and gradient or flux (Neumann). Rewriting it as $\frac{\partial u}{\partial x} = 2-u$ one can see that it can be formulated as a Neumann boundary coefficient ($\frac{\partial u}{\partial x} = g$). The second boundary condition at $x=1$ is of the usual Dirichlet prescribed value type.
This example problem can be set up and implemented in the FEATool
MATLAB m-script language as follows. First, a 1D domain is defined with
x
as the space coordinate name. The dependent variable is
named u
and assigned first order linear FEM shape
functions sflag1
. Furthermore, a
one dimensional mesh with 10 grid cells spanning $x=0$ to $x=1$ is
created with the linegrid
command.
fea.sdim = {'x'};
fea.dvar = {'u'};
fea.sfun = {'sflag1'};
fea.grid = linegrid(10,0,1);
Although one could use the pre-defined convection and diffusion,
Poisson, or custom equation modes to define the governing PDE
equation, here it is manually and explicitly defined as ux_x +
uux_t - u_t = exp(2x)
. The
parseeqn
function checks and
parses the equation to the weak finite element syntax stored in the
fea.eqn field. For convenience a coefficient expression is also
defined for the analytic reference solution.
fea.eqn = parseeqn( 'ux_x + u*ux_t - u_t = exp(2*x)', ...
fea.dvar, fea.sdim );
fea.coef = {'u_ref', 'exp(x)'};
The _x
and _t
underscore syntax means that
the terms will be evaluated and assembled implicitly as bilinear
forms. If the underscores (and _t
test functions) are
omitted the terms will be evaluated explicitly and the problem will
become more non-linear.
The Dirichlet boundary condition should be defined for $x=1$, the second and last boundary point, where the exact solution is prescribed (in the fea.bdr.d boundary coefficient field, the first entry is left empty as no Dirichlet conditions is to be applied there).
fea.bdr.d = {[],'u_ref'};
As the fist point at $x=0$ the Robin boundary 2-u
condition is defined as a Neumann condition (in the fea.bdr.n
boundary coefficient field).
fea.bdr.n = {'2-u',[]};
Lastly, parsing, solving and postprocessing can be done with the commands
fea = parseprob(fea);
fea.sol.u = solvestat(fea);
h1 = postplot( fea, 'surfexpr', 'u', 'color', 'b' );
h2 = postplot( fea, 'surfexpr', 'u_ref', ...
'linestyle', '--', 'color', 'r' );
legend( [h1(1),h2(1)], 'computed solution', ...
'analytic solution' )
xlabel('x')
ylabel('u')
grid on
The plot and image below clearly shows how the computed solution aligns with the analytic reference solution.
Note that the solver will call the
bdrneu
FEATool function during
the solution process to analyze the fea.bdr.n field, and try to
automatically split Robin conditions into Neumann and Dirichlet parts
for separate treatment. In this example the Neumann contribution
2
will be assembled and added explicitly to the right
hand side, while the Dirichlet part -u
will added
implicitly to the corresponding boundary nodes of the stiffness matrix.
The limitation for FEATool to be able to parse Robin conditions is that the linear Dirichlet part must be easily identifiable. If for example the Robin expression is implemented as a coefficient expression, then a Robin dependent variable cannot be identified and the whole expression will be evaluated explicitly. This can be seen by making the following changes to the model problem, where although an identical solution is achieved, the solver takes more iterations to converge as the problem is more nonlinear with the explicit evaluation of the Robin boundary condition.
fea.bdr.n = {'robin_bc',[]};
fea.coef = {'u_ref', 'exp(x)';
'robin_bc', '2-u'};
fea.sol.u = solvestat( fea, 'maxnit', 50 );
postplot( fea, 'surfexpr', 'u', 'color', 'g', 'linestyle', '-.' );
This model example is available as ex_robinbc1.m
in the
FEATool examples directory
Although the model example here has been implemented using the FEATool functions and MATLAB m-script language on the command line, Robin boundary can similarly be used in the FEATool GUI by entering valid Robin expressions in the corresponding Neumann or flux boundary coefficient edit fields.