Page 331 - Process Modelling and Simulation With Finite Element Methods

318 Process Modelling and Simulation with Finite Element Methods

Boundary 24 is the “source” boundary. Since the lm computation gives it a
negative boundary integral, we should interpret it as the flux out of the domain
across that boundary. The other two methods oppose the sign of the lm method
in every instance, they have the interpretation of the flux into the domain across
that boundary. Slight numerical differences occur due to the (nx,ny) method and
“by hand” (7.3.2) having errors on the order of the grid scale. The final row is
the sum of all the boundary fluxes and is consistent only for the lm method
giving exactly zero to five decimal places. The others have cumulative errors on
the order of the grid scale. Conservation of electric charge should give flux in
equals flux out, or net flux is zero. Clearly the percentage change upon refining
the mesh by approximately four-fold the number of elements results in an order
of magnitude less change in the Im method flux estimates than in the direct
computation. As the values of the different methods are approaching upon
refining the mesh, it is clear that the Lagrange Multiplier estimate of the flux is
substantially better than the direct computation. The Lagrange multiplier is an
“integrated balance” for the constraint, and in FEM, integrated quantities are
better approximated than differentiated quantities in general. This is a feature of
the weak formulation of the PDE.
The FEMLAB 2.3 User’s Guide and Introduction [l, p. 1-4001 gives a laundry
list of caveats for the use of weak boundary constraints. We reproduce them
here for completeness, and, on advice from COMSOL, update them, now that we
have a concrete example for discussion:
Strong and weak constraints should not be mixed on adjacent
boundaries, i.e. those sharing common nodes.
You must always have a constraint on boundaries when you enable the
weak boundary constraint. N.B. only Dirichlet-type boundary
conditions count as a constraint. Neumann conditions, being natural to
FEM, even if inhomogeneous, do not count as a constraint.
Scale your equations so that all coupled quantities are the same order,
to avoid convergence difficulties. Automatic scaling of variables, a
solver parameters option, does this by default.
Discontinuous constraints are only satisfied by theoretically infinite
Lagrange multipliers. In practice, this leads to large oscillations.
Be careful not to use different element shape types between boundary
and application modes. Derivative only boundary conditions should
have lower order elements (same shape) than the “bulk.”
Iterative solvers do not like the structure of the matrices (not sparse
enough) so use incomplete LU factorization as the preconditioner for
the iterative solver.