level set method reinitialization (水平集方法的重初始化）

    1  回顾

       水平集方法LSM使用了背景网格。不同于VOF，MPM等采用背景网格节点插值界面的构造方式，在于LSM采用 phi 函数隐式描述界面， phi函数是从微分几何过来的，图像学里面叫，曲线演化。LSM对于几何信息（诸如法向量，曲率）描述精确。界面演化由全局演化方程描述，（即上一篇中使用的演化方程，并不局限解界面演化，实际上是用来更新全域演化信息。所以它需要将已知在界面上的速度场扩展到全局的速度场，当然为了计算节省。采用窄带LSM等，就是限制计算域在界面附近的区域，而不跟新远场信息）

      演化方程迭代若干步，若水平集函数梯度 phi_grad 变化剧烈，就需要reinitialization.  一般通过求解  d phi/ dt = sig(phi_0) (1 - norm(phi_grad)) 得到新 phi。该方程的稳态解自动满足 norm(phi_grad) = 1 ，即意味着在重初始化计算中，界面不用改变。这正是所希望的。

      Kevin T. Chu LSMLIB    

      UCB  LSM

   2  重初始化方程

      d  phi/ dt + sig(phi_0) ( 1- norm(phi_grad)) = 0  该方程形式类似上一篇的演化方程，故可采用Godunov格式解。需要注意的是 符号函数 sig(phi_0) 的离散化：

      if ( phi < -eps) sig(phi_0) = -1;

      elseif(phi > eps)  sig(phi_0) = 1;

      else  sig(phi_0) = 0 ;    

      重初始化方程重写如下，

       d phi / dt   +  w . phi_grad = sign(phi_0）

       w 是单位法向，由0水平集指向 phi >0 的区域。这是一个有源对流方程，可采用相关数值解法。速度w = 1。 那么， 求解至 T 时刻， DSF phi(x, T)给出了到边界距离不超过为 T的所有网格节点的水平集值。

       

      2.1 PDE重初始化方程数值解

      重初始化的基本要求，是新 phi =0 与 原 phi =0 描述同一界面线。即要求在重初始化过程中，界面保持不变。当界面经过网格点，符号函数为0，自然满足；但是当界面穿过网格线，对于一维情况将出现 sig(phi(i-1))，sig(phi(i))， 均不为0. 此时的重初始化方程中  d phi /dt 项将驱动 界面更新，而不能保持不变。     

      即1D 情况下， 当界面附近出现  phi(i-1) . phi(i) < 0 或者  phi(i) . phi(i+1) < 0，上文所给的迎风格式（Godunov) 不能正确应用。 或者说， 上文的迎风格式，实际上只能更新界面点位于网格点这种情况。 对于更普遍的情况， 必须采用单侧差分。

      此时 backward 差分形式为  phi_x ^ -   = (phi - 0) / norm (D_i)    % D 是符号距离， 表示 i/i-1 网格节点到界面的符号距离

      此时sig(phi_i) = D_i / h   % h 是网格间距

      整理方程如下：   phi_i^(n+1) = phi_i^n - dt / h { (sig(D_i) norm(phi_i^n) - D_i }

      关于sig(D_i) 的处理，一般默认 sig(D_i) = sig(phi_i^n)

      即方程如下：

      phi_i^(n+1) = phi_i^n - dt / h { (sig(phi_i^n) norm(phi_i^n) - D_i }

      这里D_i 采用三次拟合得到， phi(i-2),  phi(i-1), phi(i), phi(i+1) 四点三次插值函数 p(x), 求解在 x(i-1） < x < x(i) 区间的解，即位界面点x_p。同时 距离 d_i =  x(i) -  x_p

 

       2.2  插值型重初始化