曲面の面積と変分
下記の詳細な導出はこちら
曲面の面積を与える関数は最も基本的な汎関数のひとつで、簡単に
(1)
\begin{align} a=\int_a{da} \end{align}
と書けます。その変分（形状が変化したときの面積の微小変化)は
(2)
\begin{align} \delta a=\int_a{\mathrm{div}\delta\boldsymbol{u}\mathrm{da}} \end{align}
と表せます。
$\delta\boldsymbol{u}$は仮想変位場、$\mathrm{div}\delta\boldsymbol{u}$はその発散を表します。
一般的な膜に対する仮想仕事の原理は次のように表現できます:
(3)
\begin{align} \delta w=\int_a{\boldsymbol{\sigma}:\delta\boldsymbol{u}\otimes\nabla da} \end{align}
また、次のようにも書きますが、同じことです:
(4)
\begin{align} \delta w=\int_a{\boldsymbol{\sigma}:\frac{\partial \delta\boldsymbol{u}}{\partial \boldsymbol{x}}da}, \end{align}
ここで、
(5)
\begin{align} \delta \boldsymbol{u}\otimes \nabla=\frac{\partial \delta\boldsymbol{u}}{\partial \boldsymbol{x}}\equiv \frac{\partial \delta\boldsymbol{u}}{\partial \theta^i}\otimes \boldsymbol{g}^i. \end{align}
はベクトルの勾配(後形)、と呼ばれる2階テンソルです。対称ではありませんので、前形と異なることに注意が必要です。
膜構造ではしばしば、標準的な応力場として一様で等方的な応力場を採用します。石鹸膜はその液体としての性質から、そのような応力場が実現していると考えられ、その形成する曲面は等張力曲面と呼ばれます。
さて、そのような応力場を考察する場合、解析的には応力テンソル$\boldsymbol{\sigma}$を
(6)
\begin{align} \boldsymbol{\sigma}=\hat{\sigma}\boldsymbol{I}\equiv\hat{\sigma}\boldsymbol{g}_1\otimes\boldsymbol{g}^1+\boldsymbol{g}_2\otimes\boldsymbol{g}^2 \end{align}
に置き換えます。$\boldsymbol{I}$は単位テンソル、$\hat{\sigma}$は膜の断面に作用する単位長さ辺りの力の大きさで、単位は$N/m$に準じます。
このとき等張力曲面に対する仮想仕事の原理は次のように表現されます。
(7)
\begin{align} \delta w=\int_a{\boldsymbol{I}:\delta\boldsymbol{u}\otimes\nabla \mathrm{da}}. \end{align}
これは、次の簡潔な形式に書き換えることができます。
(8)
\begin{align} \delta w=\int_a{\delta\boldsymbol{u}\cdot\nabla \mathrm{da}}=\int_a{\mathrm{div}\delta\boldsymbol{u}\mathrm{da}}. \end{align}
なぜならば、
(9)
\begin{align} \delta\boldsymbol{u}\cdot\nabla\equiv\mathrm{div}\delta\boldsymbol{u}\equiv\frac{\partial\delta\boldsymbol{u}}{\partial\theta^i}\cdot\boldsymbol{g}^i, \end{align}
であり、また、テンソルの基本的な計算より
(10)
\begin{align} \boldsymbol{I}:\delta\boldsymbol{u}\otimes\nabla=\boldsymbol{g}_i\otimes\boldsymbol{g}^i:\frac{\partial\delta\boldsymbol{u}}{\partial\theta^\alpha}\otimes\boldsymbol{g}^\alpha=\boldsymbol{g}^i\otimes\boldsymbol{g}_i:\frac{\partial\delta\boldsymbol{u}}{\partial\theta^\alpha}\otimes\boldsymbol{g}^\alpha \end{align}
さらに
(11)
\begin{align} =\left(\boldsymbol{g}^i\cdot\frac{\partial\delta\boldsymbol{u}}{\partial\theta^\alpha}\right)\left(\boldsymbol{g}_i\cdot\boldsymbol{g}^\alpha\right)=\left(\boldsymbol{g}^i\cdot\frac{\partial\delta\boldsymbol{u}}{\partial\theta^i}\right)=\mathrm{div}\delta\boldsymbol{u} \end{align}
だからです。
ここで、$\mathrm{div}\delta \boldsymbol{u}$はベクトルの発散と呼ばれ、スカラー値です。
さて、曲面の面積の変分が
(12)
\begin{align} \delta a=\int_a{\mathrm{div}\delta\boldsymbol{u} \mathrm{da}} \end{align}
として与えられますので、よく知られた関係
(13)
\begin{align} \delta a=0\leftrightarrow\delta w=0, \end{align}
を得ます。これはつまり等張力曲面と極小曲面が本質的に同一であることを示しています。
面積汎関数の停留条件と、等張力曲面の釣り合い式の導出は、本質的に同一です。
基本的な方針は仮想仕事の原理から$\delta\boldsymbol{u}$を一番外側に括りだし、これを外すというものです。
このような操作により、一般的に表現されていた問題が、設定されたパラメータにより陽に表示され、て計算や数値計算で扱いやすい形式へと書き換えられます。
$\delta\boldsymbol{u}$を括りだす共通の手続きは部分積分で与えられます。等張力局面の仮想仕事の原理に部分積分を適用すると
(14)
\begin{align} \delta w=\int_a{\boldsymbol{I}:\delta\boldsymbol{u}\otimes \nabla da}= -\int_a{\left(\boldsymbol{I}\cdot \nabla\right)\cdot\delta\boldsymbol{u}da}+\int_l{\left(\boldsymbol{I}\cdot\boldsymbol{n}\right)\cdot\delta\boldsymbol{u}dl} \end{align}
となります。ここで
(15)
\begin{align} \boldsymbol{I}\cdot\nabla\equiv\mathrm{div}\boldsymbol{I}\equiv \frac{\partial \boldsymbol{I}}{\partial \theta^i}\cdot\boldsymbol{g}^i \end{align}
です。単位行列が発散をもつのは奇妙に見えます。しかし、曲面上の単位テンソルは曲面上の接平面表します。接平面は曲面の曲率に従って様々に向きを変えますから単位テンソルは発散を持つのです。
2階のテンソルの発散はベクトルとなります。
つまり、3次元ユークリッド空間では
(16)
\begin{align} \mathrm{div}\left(\boldsymbol{g}_1\otimes\boldsymbol{g}^1+\boldsymbol{g}_2\otimes\boldsymbol{g}^2+\boldsymbol{g}_3\otimes\boldsymbol{g}^3\right)=0, \end{align}
であるが、一方で3次元ユークリッド空間に埋め込まれた曲面上では
(17)
\begin{align} \mathrm{div}\left(\boldsymbol{g}_1\otimes\boldsymbol{g}^1+\boldsymbol{g}_2\otimes\boldsymbol{g}^2\right)\neq 0. \end{align}
であるということです。
(14)式を内部仮想仕事と解釈し、曲面に作用する力や境界に作用する力と釣り合っていると考えると、
(18)
\begin{align} \delta W_i=\delta W_e\Leftrightarrow -\int_a{\left(\boldsymbol{I}\cdot \nabla\right)\cdot\delta\boldsymbol{u}da}+\int_l{\left(\boldsymbol{I}\cdot\boldsymbol{n}\right)\cdot\delta\boldsymbol{u}dl} = \int_a{\boldsymbol{f_s}\cdot\delta\boldsymbol{u}da}+\int_l{\boldsymbol{f}_e\cdot\delta\boldsymbol{u}dl} \end{align}
境界のみ仮想変位を与えたり、曲面のみ仮想変位を与えたりできますから、$\delta\boldsymbol{u}$が外れて、（このことを短く$\delta\boldsymbol{u}$の任意性より、と書きます）
(19)
\begin{align} -\left(\boldsymbol{I}\cdot \nabla\right)=\boldsymbol{f}_s\Leftrightarrow\left(\boldsymbol{I}\cdot \nabla\right)+\boldsymbol{f}_s=\vec{\boldsymbol{0}} \end{align}
(20)
\begin{align} \left(\boldsymbol{I}\cdot \boldsymbol{n}\right)=\boldsymbol{f}_e \end{align}
これらは等張力曲面の膜応力と物体力の力の釣り合い式、境界における応力と境界力の適合条件です。右辺は外力を表しますが、それぞれ単位面積当たりの物体力、単位長さ当たりの境界力を現します。もちろん、$\boldsymbol{\sigma}=\tilde{\sigma}\boldsymbol{I}$としたときは
(21)
\begin{align} \tilde{\sigma}\left(\boldsymbol{I}\cdot \nabla\right)+\boldsymbol{f}_s=\vec{\boldsymbol{0}} \end{align}
(22)
\begin{align} \tilde{\sigma}\left(\boldsymbol{I}\cdot \boldsymbol{n}\right)=\boldsymbol{f}_e \end{align}
となります。
さて、境界を完全に固定したとしましょう。すると(14)式第二項は消去することができます。
(23)
\begin{align} \delta w=-\int_a{\mathrm{div}\boldsymbol{I}\cdot\delta\boldsymbol{u}da}. \end{align}
負号が付与されているのは、膜応力が（引っ張りの時）面積が小さくなるときに仕事をするためです。
このようにして、境界を固定された曲面について
(24)
\begin{align} \delta w=0\leftrightarrow \delta a=0\leftrightarrow\mathrm{div}\boldsymbol{I}=\vec{0} \end{align}
が得られます。これは面積汎関数の停留条件、もしくは等張力曲面の自己釣り合い条件です。
これは、等張力曲面と極小曲面が本質的に同一で、そのような曲面は具体的には平均曲率が至る所0であるような曲面であることを表しています。次節では、$\mathrm{div}\boldsymbol{I}$が平均曲率ベクトルを表していることの証明を行います。
In this section, we prove that the divergence of a unit tensor of a surface is the double of the mean curvature vector.
When this is proved, consequently, we can confirm that a minimal surface is a surface which its mean curvature is zero everywhere.
We start with the Taylor expansion of a surface, i.e.,
(25)
\begin{align} \boldsymbol{x}(\theta^1+d\theta^1,\theta^2+d\theta^2)=\boldsymbol{x}(\theta^1,\theta^2)+\boldsymbol{I}\cdot d\boldsymbol{x}+\frac{1}{2}d\boldsymbol{x}\cdot\boldsymbol{H}\cdot d\boldsymbol{x}+\cdots. \end{align}
Here, $d\boldsymbol{x},\boldsymbol{I},\boldsymbol{H}$ are respectively defined by
(26)
\begin{align} d\boldsymbol{x}\equiv d\theta^i\boldsymbol{g}_i \end{align}
(27)
\begin{align} \boldsymbol{I}\equiv \boldsymbol{x}\otimes\nabla=\frac{\partial \boldsymbol{x}}{\partial \theta^i}\otimes\boldsymbol{g}^i=\boldsymbol{g}_i\otimes\boldsymbol{g}^i=\boldsymbol{g}^i\otimes\boldsymbol{g}_i \end{align}
(28)
\begin{align} \boldsymbol{H}\equiv \nabla\otimes\left(\boldsymbol{x}\otimes\nabla\right)=\boldsymbol{g}^\alpha\otimes\frac{\partial}{\partial \theta^\alpha}\left(\boldsymbol{g}_i\otimes\boldsymbol{g}^i\right). \end{align}
The derivation is written in Taylor expansion of a surface.
Now, an arbitrary tangent vector on a surface can be represented by $\boldsymbol{I}\cdot d\boldsymbol{x}=d\boldsymbol{x}$; hence when we omit the higher terms, we obtain
(29)
\begin{align} \boldsymbol{x}(\theta^1+d\theta^1,\theta^2+d\theta^2)\simeq\boldsymbol{x}(\theta^1,\theta^2)+d\boldsymbol{x}. \end{align}
This represents a tangent plane.
We define
(30)
\begin{align} \boldsymbol{N}\equiv\boldsymbol{g}_1\cdot\boldsymbol{H}\cdot\boldsymbol{g}^1+\boldsymbol{g}_2\cdot\boldsymbol{H}\cdot\boldsymbol{g}^2=\boldsymbol{g}_i\cdot\boldsymbol{H}\cdot\boldsymbol{g}^i. \end{align}
Let us rewrite the right-hand side as
(31)
\begin{align} \mathrm{tr}\boldsymbol{H}\equiv\boldsymbol{g}_i\cdot\boldsymbol{H}\cdot\boldsymbol{g}^i. \end{align}
Such an operation is called contraction. Particularly, the above operation is called a trace.
Interestingly, the following relation holds:
(32)
\begin{align} \mathrm{div}\boldsymbol{I}=\boldsymbol{N}. \end{align}
This can be verified by comparing
(33)
\begin{align} \mathrm{div}\boldsymbol{I}\equiv\frac{\partial \boldsymbol{I}}{\partial \theta^i}\cdot\boldsymbol{g}^i \end{align}
and
(34)
\begin{align} \boldsymbol{g}_i\cdot\boldsymbol{H}\cdot\boldsymbol{g}^i=\boldsymbol{g}_i\cdot\boldsymbol{g}^\alpha\otimes\frac{\partial\boldsymbol{I}}{\partial \theta^\alpha}\cdot\boldsymbol{g}^i=\frac{\partial \boldsymbol{I}}{\partial \theta^i}\cdot\boldsymbol{g}^i. \end{align}
In summary, the following relations holds:
(35)
\begin{align} \boldsymbol{N}=\mathrm{div}\boldsymbol{I}=\mathrm{tr}\boldsymbol{H}. \end{align}
It is very interesting that $\boldsymbol{N}$ represents a normal (orthogonal to the tangent plane) vector and its norm is the double of a quantity called mean curvature.
We first prove that it represents a normal vector.
The calculation is like as follows:
(36)
\begin{align} \boldsymbol{N}=\frac{\partial \boldsymbol{I}}{\partial \theta^i}\cdot\boldsymbol{g}^i=\frac{\partial \left(\boldsymbol{g}_i\otimes\boldsymbol{g}^i\right)}{\partial \theta^\alpha}\cdot\boldsymbol{g}^\alpha=\frac{\partial\left( \boldsymbol{g}^i\otimes\boldsymbol{g}_i\right)}{\partial \theta^\alpha}\cdot\boldsymbol{g}^\alpha=\frac{\partial\boldsymbol{g}^i}{\partial \theta^\alpha}\otimes\boldsymbol{g}_i\cdot\boldsymbol{g}^\alpha+\boldsymbol{g}^i\otimes\frac{\partial \boldsymbol{g}_i}{\partial \theta^\alpha}\cdot\boldsymbol{g}^\alpha \end{align}
Hence,
(37)
\begin{align} \boldsymbol{N}=\frac{\partial \boldsymbol{g}^i}{\partial \theta^i}+\boldsymbol{g}^i\otimes\frac{\partial \boldsymbol{g}_i}{\partial \theta^\alpha}\cdot\boldsymbol{g}^\alpha \end{align}
ここで、少し複雑な法則を用います。共変基底、反変基底の関係$\boldsymbol{g}_i\cdot\boldsymbol{g}^j=\delta_{i}^{\cdot j}$を思い出しましょう。いずれの場合も0か1ですから、任意の共変基底と反変基底、任意の座標パラメータ$\theta^\alpha$について
(38)
\begin{align} \frac{\partial\left( \boldsymbol{g}_i\cdot\boldsymbol{g}^j\right)}{\partial \theta^\alpha}=0 \end{align}
が成り立ちます。従って、
(39)
\begin{align} \frac{\partial \boldsymbol{g}_i}{\partial \theta^\alpha}\cdot \boldsymbol{g}^j+\boldsymbol{g}_i\cdot\frac{\partial \boldsymbol{g}^j}{\partial \theta^\alpha}=0 \end{align}
さらに、
(40)
\begin{align} \frac{\partial \boldsymbol{g}_i}{\partial \theta^\alpha}\cdot \boldsymbol{g}^j=-\boldsymbol{g}_i\cdot\frac{\partial \boldsymbol{g}^j}{\partial \theta^\alpha} \end{align}
従って、$\alpha=j$の特別な場合についても
(41)
\begin{align} \frac{\partial \boldsymbol{g}_i}{\partial \theta^j}\cdot \boldsymbol{g}^j=-\boldsymbol{g}_i\cdot\frac{\partial \boldsymbol{g}^j}{\partial \theta^j} \end{align}
が成り立ちます。代入して
(42)
\begin{align} \boldsymbol{N}=\frac{\partial \boldsymbol{g}^i}{\partial \theta^i}-\boldsymbol{g}^j\otimes\frac{\partial \boldsymbol{g}^i}{\partial \theta^i}\cdot\boldsymbol{g}_j \end{align}
さて、第二項ですが、$\frac{\partial \boldsymbol{g}^i}{\partial \theta^i}\cdot\boldsymbol{g}_j$がスカラーですから、
(43)
\begin{align} \left(\frac{\partial \boldsymbol{g}^i}{\partial \theta^i}\cdot\boldsymbol{g}_j\right)\boldsymbol{g}^j \end{align}
と書けます。従って、
(44)
\begin{align} \boldsymbol{N}=\frac{\partial \boldsymbol{g}^i}{\partial \theta^i}-\left(\frac{\partial \boldsymbol{g}^i}{\partial \theta^i}\cdot\boldsymbol{g}_j\right)\boldsymbol{g}^j \end{align}
となります。
単位法線ベクトル
(45)
\begin{align} \boldsymbol{n} \end{align}
を定義します。大きさが1で、接平面と直交していますから、共変基底、反変基底の区別を付ける必要がありません。
この単位法線ベクトルを用いると空間の単位テンソル
(46)
\begin{align} \boldsymbol{D}=\boldsymbol{g}_1\otimes\boldsymbol{g}^1+\boldsymbol{g}_2\otimes\boldsymbol{g}^2+\boldsymbol{n}\otimes\boldsymbol{n} \end{align}
を記述できます。
曲面の単位テンソル
(47)
\begin{align} \boldsymbol{I}=\boldsymbol{g}_1\otimes\boldsymbol{g}^1+\boldsymbol{g}_2\otimes\boldsymbol{g}^2 \end{align}
と比べてみると、両者の関係がよくわかります。つまり、
(48)
\begin{align} \boldsymbol{D}-\boldsymbol{I}=\boldsymbol{n}\otimes\boldsymbol{n} \end{align}
です。すると
(49)
\begin{align} \boldsymbol{N}=\frac{\partial \boldsymbol{g}^i}{\partial \theta^i}(\boldsymbol{D}-\boldsymbol{I}) \end{align}
と書けます。他に
(50)
\begin{align} \boldsymbol{N}=\frac{\partial \boldsymbol{g}^i}{\partial \theta^i}(\boldsymbol{n}\otimes\boldsymbol{n}) \end{align}
や
(51)
\begin{align} \boldsymbol{N}=\left(\frac{\partial \boldsymbol{g}^i}{\partial \theta^i}\cdot\boldsymbol{n}\right)\boldsymbol{n} \end{align}
とかけますが、いずれも同じことを表しています。つまり
(52)
\begin{align} \frac{\partial \boldsymbol{g}^1}{\partial \theta^1}+\frac{\partial \boldsymbol{g}^2}{\partial \theta^2} \end{align}
から接平面成分が削除され、法線成分のみ取り出したものが$\boldsymbol{N}$です。
以上より$\boldsymbol{N}$は法線ベクトルであることが示されました。
さてこの法線成分について
(53)
\begin{align} \frac{\partial \boldsymbol{g}^i}{\partial \theta^i}\cdot \boldsymbol{n}=-\boldsymbol{g}^i\frac{\partial\boldsymbol{n}}{\partial \theta^i}=-g^{ij}\boldsymbol{g}_j\frac{\partial\boldsymbol{n}}{\partial \theta^i} \end{align}
について計算を続けます。
第二基本形式
(54)
\begin{align} \mathrm{II}\equiv\frac{1}{2}h_{ij}d\theta^i d\theta^j \end{align}
において
(55)
\begin{align} h_{ij}=\frac{\partial \boldsymbol{g}_j}{\partial \theta^i}\cdot\boldsymbol{n}=-\boldsymbol{g}_i\cdot\frac{\partial \boldsymbol{n}}{\partial \theta^j} \end{align}
だったことを思い出すと、
(56)
\begin{equation} =g^{ij}h_{ij} \end{equation}
と書けます。
(57)
\begin{align} \frac{\partial g^{ij}\boldsymbol{g}_i}{\partial \theta^j}\cdot\boldsymbol{n}=\left(\frac{\partial g^{ij}}{\partial \theta^j}\boldsymbol{g}_i\right)\cdot\boldsymbol{n}+h_{ij}g^{ij} \end{align}
として、第一項が零であることを用いても同じ結果を得ます。
さて、曲面上の曲線と主曲率で平均曲率は
(58)
\begin{align} H=\frac{1}{2}(\lambda_1+\lambda_2)=\frac{LG-2MF+EN}{2(EG-F^2)} \end{align}
として与えられましたが、同時に
(59)
\begin{align} H=\frac{g^{ij}h_{ij}}{2} \end{align}
とも書けました。
確かに$g^{ij}h_{ij}$は平均曲率の2倍を表しています。
面積汎関数の変分は
(60)
\begin{align} \delta a=2\int_{a}{\mathrm{div}\delta\boldsymbol{u}\mathrm{da}}. \end{align}
として与えられました。
今度はこの変分問題を最急降下法で解くための離散化をGalerkin法に基づいて行います。
Galerkin法は一般に、有限要素法の基礎、として知られていますが、実は最急降下法の基礎にもなりうることはあまり知られていません。
最終的に変分問題$\delta a=0$は
(61)
\begin{align} \sum_{j}{\nabla S_{j}}=\vec{\boldsymbol{0}}, \end{align}
の形式に離散化されます。$S_j$は要素面積です。
初めに、$n$個の自由度のみ持つ 関数 $\tilde{{\boldsymbol{x}}}(x_{1},\cdots x_{n},\theta^{1},\theta^{2})$で曲面の形状を近似します。一方で、$\boldsymbol{x}(\theta^{1},\theta^{2})$の形状の自由度は無限であったことに注意してください。
このとき、仮想仕事の原理は高々$n$個の独立な仮想変位に対してしか満たすことができなくなります。
一般に$\delta\boldsymbol{u}=\delta\boldsymbol{x}$です。そして、Galerkin法では、$\delta\boldsymbol{x}$を次式で置換します。
(62)
\begin{align} \delta\tilde{{\boldsymbol{x}}}\equiv\left[\frac{\partial\tilde{\boldsymbol{x}}}{\partial x_{1}},\cdots,\frac{\partial\tilde{\boldsymbol{x}}}{\partial x_{n}},\right]\cdot\left[\begin{array}{c} \delta x_{1}\\ \vdots\\ \delta x_{n}\end{array}\right], \end{align}
このようにして、確かに$n$個の独立な仮想変位、すなわち$\frac{\partial\tilde{\boldsymbol{x}}}{\partial x_{1}}\cdots\frac{\partial\tilde{\boldsymbol{x}}}{\partial x_{n}}$が準備されました。
この仮想変位を代入することで離散化された仮想仕事の原理
(63)
\begin{align} \delta\tilde{a}=2\left(\int_{\mathrm{a}}{\mathrm{div}\frac{\partial\tilde{\boldsymbol{x}}}{\partial x_{1}}\mathrm{da}}\delta x_{1}+\int_{\mathrm{a}}{\mathrm{div}\frac{\partial\tilde{\boldsymbol{x}}}{\partial x_{2}}\mathrm{da}}\delta x_{2}+\cdots\right)=0, \end{align}
を得ます。もしくは、$n$本の連立方程式の形式で、
(64)
\begin{align} 2\int_{\mathrm{a}}{\mathrm{div}\frac{\partial\tilde{\boldsymbol{x}}}{\partial x_{1}}\mathrm{da}}=0,\cdots,\,2\int_{\mathrm{a}}{\mathrm{div}\frac{\partial\tilde{\boldsymbol{x}}}{\partial x_{n}}\mathrm{da}}=0, \end{align}
と書き換えることもできます。これは離散化された停留条件です。
面白いことに、(60)式をよく見ると、$\delta$は$\frac{\partial}{\partial x_i}$に置き換えることが可能ですので、面積の変分の手続きがそのまま、次の書き換えの根拠になってしまいます。
(65)
\begin{align} 2\int_{\mathrm{a}}{\mathrm{div}\frac{\partial\tilde{\boldsymbol{x}}}{\partial x_{1}}\mathrm{da}}=0\Leftrightarrow\frac{\partial}{\partial x_{1}}\int_{a}{\mathrm{da}}=0,\cdots\ \mathrm{and\ so\ on}. \end{align}
これは作用$\delta$が、もともと演算子$\frac{\partial}{\partial \epsilon}$として定義されたことによります。ここで、$\epsilon$は仮想変位に付加された独立パラメータです。
以上より$\delta a=0$の離散化が次のようにして得られました。
(66)
\begin{align} \sum_{j}{\nabla S_{j}}=\vec{\boldsymbol{0}}. \end{align}
