FEM for Advection-Diffusion Equation

移流拡散方程式の有限要素法による数値解法

Last Update:2009年05月15日

定常移流拡散方程式

移流拡散方程式とはスカラー量 $\phi$ が速度で輸送されつつ，拡散される方程式である。

物質と同じ速度で動く視点からみれば、輸送の現象は関係ないので、単純な拡散方程式と等価になる。よって移流拡散方程式は次のように書くことができる。

$\left.\frac{\partial\phi}{\partial t}\right|_{\b{X}} = \nabla_x \cdot (\nu\nabla_x\phi)$

$\left.\frac{\partial\phi}{\partial t}\right|_{\b{X}} = \left.\frac{\partial\phi}{\partial t}\right|_{\b{x}} + {\b{u}}\cdot(\nabla_x\phi) = \nabla_x \cdot (\nu\nabla_x\phi)$

ここで定常であるとういう条件つまり時間的に解が変化しないという条件から、空間時間導関数 $\left.\frac{\partial\phi}{\partial t}\right|_{\b{x}}$ が０であるとして

${\b{u}}\cdot(\nabla_x\phi) = \nabla_x \cdot (\nu\nabla_x\phi)$

を得るこれが定常移流拡散方程式である。

1次元定常移流拡散方程式

ひとまず一次元での定常移流拡散方程式を考えた後に多次元（２次元、３次元）に拡張する。

一次元の定常移流拡散のイメージとしては次の図のような状況を考えると分かりやすい。

上の図は一定速度で動くベルトコンベヤーの上流側と下流側に2つの温度が与えられているような図である。この2つの温度が与えられている間の温度の分布は定常移流拡散方程式に従うと考えられる。

一次元定常移流拡散方程式の一般解

上式を一次元化すると、一次元の移流拡散方程式は次のとおり

$u\frac{\partial \phi}{\partial x}=\nu\frac{\partial^2 \phi}{\partial x^2}$

長さがの領域で両端の値が固定されているとき、解は次のようにかける。

$\phi=a e^{\frac{u}{\nu}x}+b$

但し変数 a,b は境界条件を満たすように選ばれる。ここで、 $x \rightarrow Lt$ のようにを無次元の長さに変換すると解は最終的に次のようにかける。

$\phi=a e^{\frac{uL}{\nu}t}+b=a e^{Ret}+b$

但し、 $Re=\frac{uL}{\nu}$ はレイノルズ数である。（本当はuの絶対値を取るのが正しいと思われるがここでは簡単のためにそのままにしておいた）レイノルズ数が指数関数の肩にのっているので、レイノルズ数が大きくなるほど値が急激に変化することがわかる。

移流速度による解の様子の変化

境界での値を $\phi_0,\phi_1$ とすると解は次のようにかける。

$\phi=(\phi_1-\phi_0)\frac{e^{Ret}-1}{e^{Re}-1}+\phi_0$

$Re\simeq0$ の場合 $e^{Ret}$ をテーラー展開して、 $e^{Ret}\simeq 1+Ret$ であるから、 $\phi=(\phi_1-\phi_0)t+\phi_0$ となる。つまり解は線形。
の場合 $\phi=\phi_0 \quad ( \quad for \quad 0 \le t < 1 \quad ),\qquad\qquad \phi=\phi_1 \quad ( \quad for \quad t=1 \quad )$ となる。
の場合 $\phi=\phi_0 \quad ( \quad for \quad t=0 \quad ),\qquad\qquad\phi=\phi_1 \quad ( \quad for \quad 0< t\le1 \quad )$ となる。

$\phi_0=0,\phi_1=1$ のとき、これらを図にすると次のとおり

レイノルズ数を変えたときの１次元定常移流拡散方程式の解の変化

これを表示させるたgnuplotのコマンドは以下のとおり

set sample 100
set xrange [0:1]
plot (exp(100*x)-1)/(exp(100)-1) title 'Re=100',(exp(10*x)-1)/(exp(10)-1) title 'Re=10',
(exp(0.01*x)-1)/(exp(0.01)-1) title 'Re=0.01',
(exp(-10*x)-1)/(exp(-10)-1) title 'Re=-10',(exp(-100*x)-1)/(exp(-100)-1) title 'Re=-100'

差分法における風上差分

移流拡散方程式の解が基本的には指数関数であることが分かった。数値計算を行う場合は解の一次微分や２次微分を計算するためにさまざまな計算方法がある。一次微分を計算する際によく使われる中央差分では２次までの精度であるからこのような指数関数的な解には対応できない。そこで、移流項の部分で風上差分を用いた離散化がよく使われる．

一次風上差分法
$u\frac{\partial\phi}{\partial x}=\left\{\begin{array}\displaystyle{u\frac{\phi_i-\phi_{i-1}}{h}}\qquad (u>0)\\ \displaystyle{u\frac{\phi_{i+1}-\phi_i}{h}}\qquad (u<0) \end{array}\right.$

一次の風上差分は勾配を計算したい点x_iに対して1/2hだけ上流にある点における勾配を点x_iでの勾配とすることを意味している．

中央差分不安定で風上差分が安定な直感的な理由

定常移流拡散方程式は風下に行くにしたがって解の勾配が大きくなることがわかる。

下の図は u>0 つまり左から右へ風が吹いてる時の解の様子である。

中央差分と風上差分の比較

赤色の線の勾配は点iにおける中央差分を、青色の線の勾配は点iにおける風上差分を、緑色の線の勾配は点iにおける実際の勾配を表す

上の図からも明らかなように、風上差分法は実際より小さめに、中央差分は実際より多きめに勾配の大きさが評価されることがわかる。

レイノルズ数が小さい場合は $\phi$ は直線であるから、中央差分と風上差分は実際の勾配に等しくなる。
レイノルズ数が非常に大きい場合は $\left.\frac{d \phi}{d x}\right |_i=\frac{\phi_i-\phi_{i-1}}{h}=0$ となり風上差分は実際の解と一致する。しかし風下の固定境界条件付近ではずれた大きな値をとってしまう。２回微分の項がそれと釣り合わなければならないのでここの部分で大きく凸や凹なりそれが風上に伝わって風下の固定境界付近で解が振動する。

風上差分法と数値粘性

上式は中央差分を用いて

$u\frac{\partial\phi}{\partial x} = u\frac{\phi_{i+1}-\phi_{i-1}}{2h} - \frac{|u|h}{2} \frac{\phi_{i+1}-2\phi_i+\phi_{i-1}}{h^2}$

と変形できる．上式右辺第２項は風上差分を用いることによって人工的に導入された拡散項であるので人工拡散(Artificial Diffusion)または数値粘性と呼ぶ．また，人工拡散の係数を人工拡散係数と呼ぶ．

上式のように風上差分法による人工拡散係数 $\nu_{\infty}$ は

風上差分法による人工拡散係数
$\nu^{\infty} = \frac{\|u\|h}{2}$

であたえられる．これらを用いると人工粘性を導入した移流拡散方程式は中央差分を用いることによって次のように書ける．

$u\frac{\phi^{i+1}-\phi^{i-1}}{2h} - \frac{|u|h}{2}\frac{\phi^{i-1}-2\phi^i+\phi^{i+1}}{h^2} = \nu\frac{\phi^{i-1}-2\phi^i+\phi^{i+1}}{h^2}$

$u\frac{\phi^{i+1}-\phi^{i-1}}{2h} = (\nu+\nu^{\infty})\frac{\phi^{i-1}-2\phi^i+\phi^{i+1}}{h^2}$

このとき風上差分はもともとの移流拡散方程式を

$u\frac{\partial\phi}{\partial x} = (\nu+\nu^{\infty})\frac{\partial^2\phi}{\partial x^2}$

と書き換え中央差分したものと等しい．

人工粘性の最適化

移流項の影響が大きい場合は上式は精度よく厳密解と一致することが知られているが移流の影響が小さい場合に，実際の解と比べて数値解がなまされることが知られている．そこで移流の影響に応じて人工粘性の影響を変化させることで，数値解を厳密解と一致させることができる．

具体的には人工粘性 $\nu_{opt}$ を以下のように定めればよい

$\nu_{opt} = \nu_{\infty} L(\lambda), \qquad \lambda=\frac{\nu_{\infty}}{\nu_{opt}} = \frac{1}{2}\frac{|u|h}{\nu}=\frac{1}{2}Re_c$

ここで Re_c はセルレイノルズ数であり，

$L(x) = coth(x)-\frac{1}{x}$

はランジュバン関数である．

$\lambda \rightarrow {\infty}$ で $L(\lambda)=1$ であるから，このとき $\nu_{opt}=\nu_{\infty}$ となる．すなわち，移流項の大きな問題に対して１次風上差分法の解は厳密解に一致する．また， $\lambda=0$ で L(0)=0 であるから移流項の影響が小さい場合は人工粘性の値も小さくなることが分かる．しかし，たとえ移流項の影響が小さくても人工粘性を加えた解のほうが厳密解に近くなることは注意が必要である．

以上のことから，

$u\frac{\partial\phi}{\partial x} = (\nu+\nu_{opt})\frac{\partial^2\phi}{\partial x^2}$

を中央差分で離散化すると厳密解を得ることができる．

風上差分の有限要素法への応用

Bubnov-Galerkin法と中央差分法

１次元移流拡散問題にたいして，Bubnov-Galerkin法による有限要素法と中央差分による差分法が等しいことを示す．Bubnov-Galerkin法とは補間関数の形状関数と重み関数が等しいことをいう．

微分方程式

$u\frac{\partial\phi}{\partial x} = \nu \frac{\partial^2 \phi}{\partial x^2}$

${\int^1_0}{\delta\phi}u\frac{\partial\phi}{\partial x}dx = {\int^1_0}{\delta\phi}\nu \frac{\partial^2 \phi}{\partial x^2}dx$

${\int^1_0}u{\delta\phi}\frac{\partial\phi}{\partial x}dx = -{\int^1_0}\nu\frac{\partial\delta\phi}{\partial x}\frac{\partial\phi}{\partial x}dx$

Bubnov-Galerkin法により形状関数と重み関数が互いに等しい関数であらわされているとする．つまり，

$\phi = {N^i}{\phi^i}$ 、 $\delta\phi = {N^j}{\delta\phi^j}$

のように補間されているとする．上式を上式に代入すると

${\int^1_0}u{N^i\delta\phi^i}\frac{\partial N^j\phi^j}{\partial x}dx = -{\int^1_0}\nu\frac{\partial N^i\delta\phi^i}{\partial x}\frac{\partial N^j\phi^j}{\partial x}dx$

両辺から節点iでの重みの値 $\delta\phi_i$ を削除すると

${\int^1_0}u{N^i}\frac{\partial N^j\phi^j}{\partial x}dx = -{\int^1_0}\nu\frac{\partial N^i}{\partial x}\frac{\partial N^j\phi^j}{\partial x}dx$

重み関数 N^i が０でない値を持つのは x=[x_i-h,x_i+h] の範囲である．この範囲で形状関数 N^j が０でない値を持つのは j=i-1,i,i+1 の時のみである．よって重み関数 N^i についての以下の式が成り立つ．

$u{\int^{x_i+h}_{x_i-h}} {N^i}\frac{\partial N^{i-1}\phi^{i-1}}{\partial x} + {N^i}\frac{\partial N^i\phi^i}{\partial x} + {N^i}\frac{\partial N^{i+1}\phi^{i+1}}{\partial x}dx\\ = -\nu{\int^{x_i+h}_{x_i-h}} \frac{\partial N^i}{\partial x}\frac{\partial N^{i-1}\phi^{i-1}}{\partial x} + \frac{\partial N^i}{\partial x}\frac{\partial N^i\phi^i}{\partial x}dx + \frac{\partial N^i}{\partial x}\frac{\partial N^{i+1}\phi^{i+1}}{\partial x}dx$

この積分を実行すると以下の式が得られる．

$u\frac{\phi^{i+1}-\phi^{i-1}}{2h} = \nu\frac{\phi^{i-1}-2\phi^i+\phi^{i+1}}{h^2}$

これは中央差分法を用いてこの問題を離散化した式と同一である．よってBubnov-Galerkin法と中央差分法は１次元の定常移流拡散問題では同一になることが分かる．

風上差分の有限要素法への適応

中央差分の場合は適切な人工粘性 $\nu_{opt}$ を導入することで数値解を厳密解にすることができることがわかった．そこで有限要素法の場合でも人工粘性 $\nu_{opt}$ を付加することによって数値解を厳密解にすることができると考えられる．

人工粘性を付加した弱形式は次のとおり

$u\frac{\partial\phi}{\partial x} = (\nu+\nu_{opt})\frac{\partial^2 \phi}{\partial x^2}$

上式を参考にしてこの方程式の有限要素法による離散化は

${\int^1_0}{N^i}u\frac{\partial N^j\phi^j}{\partial x}dx = -{\int^1_0}(\nu+\nu_{opt}) \frac{\partial N^i}{\partial x}\frac{\partial N^j\phi^j}{\partial x}dx$

となる．この式を変形して人工粘性の項を左辺の移流項に移動させると

${\int^1_0}(N^i+\frac{\nu_{opt}}{u}\frac{\partial N^i}{\partial x})u \frac{\partial N^j\phi^j}{\partial x}dx = -{\int^1_0}\nu\frac{\partial N^i}{\partial x}\frac{\partial N^j\phi^j}{\partial x}dx$

となる．この式の意味することは，上式の左辺つまり弱形式の移流項に対して重み関数を

$N^i \quad \rightarrow \quad \bar{N}^i\equiv N^i+\frac{\nu_{opt}}{u}\frac{\partial N^i}{\partial x}$

のように変形させたことに等しい．この形状関数 $\bar{N}^i$ はもともとの形状関数 N^i に比べて風上側の重みが大きく，風下側の重みが小さくなっており，風上差分と同じ効果をもっていることが分かる．

多次元移流拡散問題の安定化

多次元における定常移流拡散方程式は次のようにかける．

${\b{u}}\cdot(\nabla_x\phi) = \nabla_x \cdot (\nu\nabla_x\phi)$

定常１次元移流拡散方程式の場合は粘性項に数値粘性係数を付加し，Bubnov-Galerkin法による離散化を行うことで解を精度よく求めることができた．しかし上式の粘性項に数値粘性係数をそのまま付加した式

${\b{u}}\cdot(\nabla_x\phi) = \nabla_x \cdot \{(\nu+\nu_{opt})\nabla_x\phi\}$

をBubnov-Galerkin法によって離散化するとなまった非現実的な解を与えることが知られている．風上差分のもともとの考えは風上方向に移流項の計算点を移動させることで流線方向に数値拡散を加えるというものだった．しかし，上式は粘性をただ大きくしただけで，流線方向に数値拡散を加えるだけではなく，流線方向に垂直な方向に対しても数値拡散を加えているので実際よりも粘性が大きくなってしまうので問題である．

そこで流線方向にのみ数値拡散項を付加することを考える．そのための準備として次のような流線方向 $\b{e}_s$ の空間微分オペレータを考える．

$\frac{\partial}{\partial x_s} = {\b{e}_s} \cdot \nabla_x \otimes = \frac{\b{v}}{|{\b{v}}|} \cdot \nabla_x \otimes = \frac{v_i}{|{\b{v}}|} \frac{\partial}{\partial x_i}$

このオペレータを用いると例えば，

$\frac{\partial \mathscr{A}}{\partial x_s}$

は流線方向 $\b{e}_s$ に沿った物理量 $\mathscr{A}$ の単位長さあたりの変化量であり，流線方向の勾配(gradient)である．

上式であらわされる多次元の定常移流拡散方程式に流線方向の人工粘性係数を付け加えた式は次のように書ける

${\b{u}}\cdot(\nabla_x\phi) = \nu \nabla^2\phi + \frac{\partial}{\partial x_s} \left( \nu_{opt} \frac{\partial\phi}{\partial x_s}\right)$

${\b{u}}\cdot(\nabla_x\phi) = \nu \nabla^2 \phi + {\b{e}_s} \cdot \[ \nabla_x \{\nu_{opt}{\b{e}_s}\cdot (\nabla_x\phi)\}\]$

上式を弱形式化すると

$\int_v {\delta\phi}{\b{u}}\cdot(\nabla_x\phi) dv = -\int_v \nu \nabla_x{\delta\phi} \cdot \nabla_x\phi + {\delta\phi}{\b{e}_s} \cdot \[ \nabla_x\{\nu_{opt}{\b{e}_s}\cdot(\nabla_x\phi)\}\] dv$

右辺第３項について変形を行う．

$\int_v \delta\phi{\b{e}_s} \cdot \[ \nabla_x\{\nu_{opt}{\b{e}_s}\cdot(\nabla_x\phi)\}\] dv\\ = \int_v {\b{e}_s} \cdot \[ \nabla_x \{ \nu_{opt}{\delta\phi}{\b{e}_s}\cdot(\nabla_x\phi)\} \] - \{{\b{e}_s}\cdot(\nabla_x{\delta\phi})\} \{\nu_{opt}{\b{e}_s}\cdot(\nabla_x\phi)\} dv\\ = \int_v \nabla_x \cdot [ \{{\delta\phi}\nu_{opt}{\b{e}_s}\cdot(\nabla_x\phi)\} {\b{e}_s}\] - \underbrace{\left\{\nabla_x \cdot {\b{e}_s}\right\}}_{=0} \{{\delta\phi}\nu_{opt}{\b{e}_s}\cdot(\nabla_x\phi)\} - \{{\b{e}_s}\cdot(\nabla_x{\delta\phi})\} \{\nu_{opt}{\b{e}_s}\cdot(\nabla_x\phi)\} dv$

第１項についてガウスの発散定理を適応する．また，流線方向 $\b{e}_s$ は一定とすると，第２項は０となる．よって上式は

$\int_s{\b{n}}\cdot\underbrace{\[ \{{\delta\phi}\nu_{opt}{\b{e}_s}\cdot(\nabla_x\phi)\}{\b{e}_s} \]}_{=0}ds - \int_v\{{\b{e}_s}\cdot(\nabla_x{\delta\phi})\}\{{\nu_{opt}}{\b{e}_s}\cdot(\nabla_x\phi)\} dv\\ = -\int_v\{{\b{e}_s}\cdot(\nabla_x{\delta\phi})\} \{{\nu_{opt}}{\b{e}_s}\cdot(\nabla_x\phi)\} dv\\ = -\int_v\{\frac{\b{v}}{|\b{v}|}\cdot(\nabla_x{\delta\phi})\} \{{\nu_{opt}}\frac{\b{v}}{|\b{v}|}\cdot(\nabla_x\phi)\} dv\\ = -\int_v\{\frac{\nu_{opt}}{|\b{v}|^2}{\b{v}}\cdot(\nabla_x{\delta\phi})\} \{{\b{v}}\cdot(\nabla_x\phi)\} dv$

と式変形される．これを上式に代入すると

$\int_v \left\{{\delta\phi}+\frac{\nu_{opt}}{|\b{v}|^2}{\b{v}}\cdot(\nabla_x{\delta\phi})\right\}\left\{{\b{v}}\cdot(\nabla_x\phi)\right\} dv = -\int_v\nu \nabla{\delta\phi} \cdot \nabla{\phi}$

これは移流項の重み関数 $\delta\phi$ を $\delta\bar{\phi}$ と変化させることに等しい．

$\delta\phi \quad \rightarrow \quad \delta\bar{\phi} = {\delta\phi}+\tau{\b{v}}\cdot(\nabla_x{\delta\phi})$

但し， $\tau$ は以下で定義される安定化パラメータである．

安定化パラメータ
$\tau = \frac{\nu_{opt}}{\|\b{v}\|^2}$

ここで付加した数値拡散について詳しく調べる．数値拡散の項を変形すると，

${\b{e}_s} \cdot \left\{ \nabla( \nu_{opt} {\b{e}_s} \cdot {\nabla\phi} ) \right\} \\ = \nabla\cdot \left\{{\b{e}_s}(\nu_{opt}{\b{e}_s}\cdot{\nabla\phi})\right\} - \underbrace{(\nabla\cdot{\b{e}_s})}_{=0} (\nu_{opt}{\b{e}_s}\cdot{\nabla\phi})\\ =\nabla\cdot(\underbrace{(\nu_{opt}{\b{e}_s}\otimes{\b{e}_s})}_{Tensor \quad Viscosity}\cdot{\nabla\phi})$

流線の方向の変化は無視できるとした．ここで

$\nu_{opt}{\b{e}_s}\otimes{\b{e}_s} = \nu_{opt}\frac{\b{v}}{|\b{v}|}\otimes\frac{\b{v}}{|{\b{v}}|} = \frac{\nu_{opt}}{|\b{v}|^2}{\b{v}}\otimes{\b{v}} = \tau{\b{v}}\otimes{\b{v}}$

はテンソル粘性(Tensor Viscosity)という．流線の方向にのみ粘性となる値が存在するような異方性を有していることがわかる．

多次元移流拡散問題の有限要素法離散化

多次元問題の安定化パラメータ

安定化パラメータは次のように表されるのであった。

$\tau = \frac{\nu_{opt}}{||\b{v}||^2}$

$\nu_{opt}$ は次のようにすると、一次元において節点における解が厳密解と一致するのであった。

$\nu_{opt} = \frac{||\b{v}||h}{2}(\coth(\lambda)-\frac{1}{\lambda})$ 、 $\lambda=\frac{1}{2}\frac{||\b{v}||h}{\nu}$

ここではメッシュサイズである。1次元の場合はメッシュサイズは線要素の長さと等しいので特に問題ないのであるが、２次元３次元の場合はメッシュサイズは簡単には決まらない。２次元や３次元の要素ではは流速の方向についての要素の長さを代表していなければならない。

Tezduyarは1991年の論文*1の中で次の様に流速方向のメッシュサイズを決定している。

$h = 2 \left( \sum_{a=1}^{n} \left|\b{s} \cdot \nabla_x \b{N}_a \right|\right)^{-1}$

ここでは要素の内挿関数の数である。また $\b{s}$ は流速ベクトル $\b{v}$ を長さが１に規格化した方向ベクトルである。

また、この論文では次のような近似を用いている。

$\coth(\lambda)-\frac{1}{\lambda}\simeq z(\lambda)=\{\quad\begin{array}\frac{1}{3}\lambda\\1\end{array}\qquad\qquad\begin{array}(\quad 0\le\lambda\le3 \quad )\\ (\quad 3\lt \lambda \quad )\end{array}$

近似関数ｚと近似される前の関数は次のグラフのとおり

サンプルプログラム

問題設定

次のサンプルプログラムは大きさ１×１の矩形領域における次のような移流拡散方程式(Advection Diffusion Problem)を解く問題である。

$\b{u}\cdot\nabla\phi-0.001{\nabla}^2\phi=0 \qquad ( \quad in \quad \Omega = \{x,y:-0.5<x,y<0.5 \} \quad )$

$\b{u}=(\begin{array}x\\-y\end{array})$

$\{\begin{array} \phi=0 \qquad ( \quad on \quad \partial\Omega \quad ) \\ \phi=(\sin2\pi x)^5 \qquad ( \quad on \quad \{x,y:0.0<x<0.5,y=0.0 \} \quad) \end{array}$

矩形領域に原点を通る切り込みが入っており、そこに固定境界条件がはいっている。流速は原点を中心に円を描くように反時計周りで流れているものとする。

問題の離散化は安定化有限要素法を用いている。

ダウンロード

DOWNLOAD	GMRes.zip

OS:WindowsXP, Windows2000
開発環境:VC++ 2005
言語:C++

結果の表示

プログラムを実行すると'solution.dat'という結果ファイルが生成される。

'solution.dat'はgnuplot用のデータ形式をしている

windows版のgnuplotの導入は以下のページが詳しい

初歩 gnuplot入門: http://auemath.aichi-edu.ac.jp/~khotta/ghost/gnuplot.html

物理の鍵しっぽ(Windows版のGnuplotの初期設定): http://hooktail.org/computer/index.php?Windows%C8%C7Gnuplot%A4%CE%BD%E9%B4%FC%C0%DF%C4%EA

値を高さで表した表面を表示

gnuplotを起動した後、プログラムを実行した'solution.dat'のあるディレクトリに移動し、次のようなコマンドを打つと結果が表示される。

gnuplot> set style data line
gnuplot> splot 'solution.dat'

等高線を表示

また次のようなコマンドを打つとgnuplotは下のような等高線を表示する。

gnuplot> set nosurface
gnuplot> set contour base
gnuplot> set view 180,90
gnuplot> splot 'solution.dat'

参考文献

有限要素法による流れのシミュレーション	日本数値流体力学学会有限要素法研究委員会 (著)
流れの有限要素法解析(2)	棚橋隆彦 (著), 伊理正夫 (著), 三好俊郎 (著), 小国力 (著), 小柳義夫 (著)

*1Tezduyar, T. E. : “Stabilized finite element formulation for incompressible flow computations”, Advance in Applied Mechanics, Vol. 28 , pp. 1-44 (1991)