0

To give some background on the problem, I'm trying to solve the following PDE via finite differences:

$$\frac{\partial\sigma}{\partial t}=\frac{\partial}{\partial x}\left[ \kappa(x,t)\left(\frac{\partial\sigma}{\partial x}+G\right)\right]$$

After expanding the equation it takes the form of a convection diffusion equation:

$$\frac{\partial\sigma}{\partial t}=\kappa(x,t)\frac{\partial^2\sigma}{\partial x^2}+\frac{\partial\kappa(x,t)}{\partial x}\frac{\partial\sigma}{\partial x}+G\frac{\partial\kappa(x,t)}{\partial x}$$

It should also be noted that $k(x,t)\geq 0$, $G>0$, and $\frac{\partial\kappa(x,t)}{\partial x}\leq0$. All the coefficients vary in time and space, and in some regions the PDE is convection-dominated, which is known to make the problem more difficult. Right now, however, I've managed to get results which make physical sense using an upwind semi-discrete scheme; essentially, I discretize the equation in space using central differences on the diffusive (second order) term and forward differencing on the convective (first order) term. I then use a stiff ODE solver in MATLAB to solve the resulting system of ODEs.

The equation is also subject to Neumann conditions $\sigma_x(0,t)=\sigma_x(L,t)=-G$, with $G>0$, and this is where I'm running into some trouble. I'm using central differences to incorporate the Neumann conditions into the scheme; if $x_1,x_2,\dots,x_n$ are the mesh nodes, I use $\sigma_x(0,t)\approx (\sigma(x_2,t)-\sigma(x_0,t))/2\Delta x$ and plug the expression for $\sigma(x_0,t)$ into the scheme centered at $x_1$. A similar approach is taken for the other boundary condition.

The solution behavior seems to be somewhat accurate on the interior nodes, but at the $x=0$ boundary, the function $\sigma(0,t)$ either oscillates, or increases monotonically when it should decrease monotonically, which is a big problem. It doesn't grow without bound, but the oscillations/monotonic increases make the solution unphysical.

I have two questions: (1) Upwinding is one approach used to limit spurious oscillations at the interior nodes, but are there any techniques to limit spurious oscillations at boundaries with Neumann conditions? (2) Is there a way to impose some kind of monotonicity or non-positivity restriction on the evolution of the solution at $x=0$?

kieransquared
  • 2,750
  • What's the sign of $-G b(0,t) + F(0,t)$? (I am wondering if maybe your expectation that $\sigma(0,t)$ is a decreasing function is unrealistic.) – Ian Feb 12 '21 at 21:40
  • By the way, how fast does $a$ vary with $x$? If it varies too fast then the neglect of the $a_x$ term that arises from the divergence-form equation is not appropriate. – Ian Feb 12 '21 at 21:43
  • $F(x,t)$ is actually equal to $Gb(x,t)$ everywhere, and $a$ doesn't vary too quickly with $x$, it mainly varies with time and retains the same qualitative $x$-dependence. I left the equation in its most general form but there is some dependence between the coefficients, I'll update the question. – kieransquared Feb 12 '21 at 21:47
  • OK. Numerically, you are actually better off using the divergence form version throughout rather than expanding, especially in the convection dominated regime. You also might consider shifting everything so that you are solving for $\sigma+Gx$ (then shift back if you need to). This last remark probably won't do all that much though. – Ian Feb 12 '21 at 21:58
  • 1
    Try looking at https://math.stackexchange.com/questions/1509291/numerical-solution-of-non-constant-coefficient-diffusion-equation-via-finite-dif/1509321?noredirect=1#comment8306556_1509321 applied to $c=\sigma + Gx$ and with homogeneous Neumann boundary conditions, see if it works reasonably well. – Ian Feb 12 '21 at 22:20
  • @Ian I tried using the spatially symmetric conservative scheme that you mentioned in that other post (along with shifting to "homogenize" the BCs), and it seems to have removed some key aspects of the solution. Previously, the upwind scheme applied to the expanded form had a solution similar to a traveling wave with positive and negative components (which is what I expect) but now my solution is entirely positive and there's no wave-like behavior. Any suggestions? – kieransquared Feb 19 '21 at 15:00
  • Did you remember to backshift by $Gx$? And did you make sure to use the different values of $\kappa$ (as this is what drives the convection in this form)? If that didn't help, then I would be curious to have some idea of what your $\kappa$ actually looks like so I can play with this myself. – Ian Feb 19 '21 at 16:22
  • Yup, I backshifted by $Gx$ and made sure the initial conditions were shifted too. $\kappa$ is a nonlinear function of the solution data of another diffusion process (basically, this models stress generation due mass transport caused by electron momentum transfer in an Al(Cu) metal line), but it's sigmoidal in space and closely resembles the function $$u(x,t)=\begin{cases} 1 & x \leq f(t) \ 0 & x> f(t) \end{cases}$$ except with a smooth transition around x=f(t), and where $f(t)$ roughly grows like $\sqrt{x}$. – kieransquared Feb 20 '21 at 03:39
  • Physically, $\kappa$ depends on the amount of $Cu$ solute at a particular point, and $Cu$ impedes stress generation so it reduces the diffusivity. But as Cu gets depleted from the $x=0$ end of the line, the diffusivity $\kappa$ increases, and the Cu depletion also has a wave-like shape as its diffusion process also has a convective component. Hope this sheds some light on the problem. – kieransquared Feb 20 '21 at 03:44
  • Interesting setup. I may tinker with some numerics of my own with this. Can you explain what you meant about the growth of $f(t)$? – Ian Feb 20 '21 at 05:01
  • Also, is your mesh fine enough to resolve the interface between the low conductivity region and the high conductivity region? – Ian Feb 20 '21 at 15:39

0 Answers0