3
$\begingroup$

I am trying to solve the Lamm eqn numerically, using 1-D FEM. What I observe is spurious oscillations at one boundary ($r_b$), which quickly leads to unreasonable values over the entire range. I am looking for advice on how to modify the FE formulation to eliminate these oscillations.

The Lamm eqn is $${\partial C \over \partial t} + {1\over r}{\partial (r J) \over \partial r}=Q, \ \ r_a\leq r\leq r_b,$$ and $$J=s \omega^2 r C - D {\partial C \over \partial r}$$ where $s$ is the sedimentation coefficient, $D$ is the diffusion coefficient, and $\omega$ the angular velocity. The boundary conditions are that the fluxes $J(r_a,t) = J(r_b, t)= 0$. For simple cases, $s$ and $D$ are constants, but they may also be functions of $C$. $Q$ is a "source function", but it is 0 for our case. The initial condition is typically $C$=constant.

The current FEM formulation uses piecewise linear ("hat") functions for discretization in space, and an implicit method to compute $C(t+\Delta t)$ from $C(t)$.

EDIT: Typical values are $r_a$=5.9, $r_b$=7.2, $C_0$=1, $s=4\times 10^{-13}$, $\omega^2=2.7\times 10^7$, $D=3.6\times 10^{-7}$.

At long times, the solution is proportional to $\exp(r^2)$.

$\endgroup$
  • $\begingroup$ Could you add some typical numerical values for the different parameters? I'm guessing that $D$ is small compared with $s\omega^2$. There are various ways to deal with these oscillations. In 1D, the most basic approaches more or less boil down to making $D$ larger, somewhat arbitrarily, until the oscillations disappear. $\endgroup$ – Bill Greene Apr 10 '17 at 21:01
  • $\begingroup$ Your equation reminds me of a convection-diffusion equation. If the convection is large w.r.t to diffusion term ($D$ is relatively small), then your solution behaves like a solution of convection equation. The simplest FE methods do not work well for convection equation because for convection-dominated equations monotonicity is essential. Non-monotonic methods produces oscillations, which you might observe. So the suggestion is to try a simple upwind finite difference method, for example. $\endgroup$ – VorKir Apr 11 '17 at 1:05
  • $\begingroup$ I have added typical numerical values (cgs units). For these $s\omega^2=1\times 10^{-5}$, so it is larger than $D$. $\endgroup$ – Woody20 Apr 11 '17 at 2:49
3
$\begingroup$

The Lamm equation is a specific example of a convection-diffusion PDE. You can find many references on the difficulties of numerical solution of this equation (e.g. Leveque).

Oscillations in the solution near the boundaries, due to characteristics of the numerical method are a common problem. A key parameter in understanding the likelyhood of oscillations is the mesh Peclet number given by

$$ P = \frac{s\omega^2r\Delta x}{2D} $$ where $\Delta x$ is the mesh spacing. Oscillations almost always occur when this number is $>1$. For the Lamm equation and the constants you provided, for a reasonable mesh, this number is $<1$. However, due to the boundary conditions, particularly at the right end, there is a very steep gradient in the solution. This, along with even the moderate Peclet number is responsible for the oscillations there.

The obvious solution is to refine the mesh to reduce the Peclet number. This can be done non-uniformly so that, for example, the mesh is much more refined at the right end. Often, this produces a model that too computationally-costly to solve, but that does not appear to be the case here.

Another alternative is to artificially increase the value of $D$. One of the comments suggested to instead use a forward-difference (upwind) finite difference method. This has the unfortunate side-effect of significantly increasing $D$ so that the solution is very "smeared" compared with exact. In the finite element method, you can achieve the same effect by slightly increasing $D$ directly.

By doing a simple web search, I was able to find a number of papers related to numerical solution of the Lamm equation for the particular application you are interested in. Many of them describe the difficulties of using basic FE methods for this problem and propose more sophisticated alternatives, such a moving-grid methods, to overcome these problems. I think this one, by Liu and Chen, does a particularly nice job in this regard.

Finally, using the parameters you provided, I computed the following solution to your equation using $500$ two-node finite elements and an implicit time integration algorithm. As you can see, for these parameters and this level of mesh refinement, there is no noticeable oscillation at the right end. enter image description here

$\endgroup$
  • $\begingroup$ Tnx for looking at this. I ran the same FEM solution as you, and I agree there are no oscillations at high $r$, although there is a small oscillation at the very start, $t=0.01$. The problems I'm having trouble with have $D$ dependent on $C$ (increasing $C$ increases $D$). When I do get oscillations, decreasing the mesh spacing, say from 1500 pts to 24000 pts, doesn't seem to help. $\endgroup$ – Woody20 Apr 11 '17 at 19:41
  • $\begingroup$ Exactly how does $D$ vary with $C$? $\endgroup$ – Bill Greene Apr 11 '17 at 20:06
  • $\begingroup$ Roughly speaking, $D=D_0(1+2bC)$, where $b$ is a constant coeff. It's sometimes more complicated. $\endgroup$ – Woody20 Apr 11 '17 at 22:26
  • $\begingroup$ I tried this expression for $D$ with $b=.1, b=1, b=10$ but was unable to reproduce the problem you describe. But there is a much broader question that needs to be answered. Finite element solution of the Lamm equation for this particular application seems to be well-studied (e.g. the Liu paper I reference and many papers by Schuck). What is about their proposed approaches that is not applicable for your case? $\endgroup$ – Bill Greene Apr 12 '17 at 16:51
  • $\begingroup$ We have been using the FE method originally described in their ref [3] for many years (pre Liu and Chen), without much difficulty. Now, in a more complex situation (2 coupled C's), we have oscillations that do not respond to smaller mesh size. I would consider using the L-C method, but it will be a big undertaking to figure out the details from their brief description. Before undertaking that, I'd like some assurance that it would work. On another topic, may I ask, are you using a FE package to do your tests? $\endgroup$ – Woody20 Apr 12 '17 at 18:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.