趣味で計算流砂水理

趣味で計算流砂水理 Computational Sediment Hydraulics for Fun Learning

数値計算とか河川工学とかプログラミングのことを書いています

MENU

pythonで地域メッシュごとのポリゴンデータからGeoTIFFを作成する：国土数値情報平年メッシュデータを例に

R

TL;DR
  • 地域メッシュポリゴンデータからGeoTIFFを作る手順を解説
  • pythonのrioxarrayを使うと簡単
  • テストケースで作成した平均年間雨量のデータが興味深い。

結果だけ参照されたい方は「 国土数値情報：平均年間雨量メッシュデータ」まで読み飛ばして下さい。

モチベーション

国土数値情報 | 平年値メッシュデータ が、1次地域メッシュごとのshapefileでポリゴンデータというかなり扱いにくい形式だったので、この手のデータの近年の主流であるGeoTIFFに変換してみました。

地域メッシュとは？

過去記事 欠損したグリッド標高テキストデータからgeotiffを作成するpythonスクリプト - 趣味で計算流砂水理 Computational Sediment Hydraulics for Fun Learning で解説した航測測量データの図郭は平面直角座標系（XY）で定義されますが、地域メッシュは地理座標系（緯度経度）で下表のとおりに定義されます。（出典はこちら） 左下を起点に番号が振られています。

参照：統計局ホームページ/地域メッシュ統計 の地域メッシュ統計の特質・沿革

元shapefileファイルのデータ構造

shapefileファイルのデータ構造の主な特徴は次のとおりです。（詳細はこちらを参照）

  • 1次メッシュごとに153個のファイル作成されており、各ファイルに3次メッシュごとのポリゴンが格納されています。
  • そのため、各ファイルのポリゴン数は80X80になりますが、水域はデータが作成されていないため、ファイルごとにポリゴン数が異なります。
  • 各ポリゴンの属性に3次メッシュコードが格納されています。

GeoTIFF用二次元配列の作成

80x80の二次元配列を作成して3次メッシュコードから値を入力します。 型はuint16（範囲は0～65535）を使用しますので、nodataを65000にしています。

こちらのG02-12_5337-jgd_GMLのファイルから平均年間雨量を雨量のデータからGeoTIFFを作ります。

import geopandas as gpd
import numpy as np

fp = "G02-12_5337-jgd.shp"
gdf = gpd.read_file(fp)

ind = gdf["G02_001"].values # 3rd mesh code
val = gdf["G02_014"].values # Total annual rainfall

# size 1st order Chiki-Mesh
x = np.full((80,80), int(65000), dtype='uint16') # y, x

for dp, v in zip(ind, val):
    # yindex
    iy = 10 * int(dp[4]) + int(dp[6])
    # yindex
    ix = 10 * int(dp[5]) + int(dp[7])
    x[iy,ix] = int(v)

メッシュセンターの緯度経度の設定

前述のとおり、メッシュコードからポリゴン左下の座標を計算できます。 値はメッシュのセンターで定義するため緯度経度方向ともに半メッシュずらします。

dp = ind[0] #3次メッシュコード
lonp = float(100) + float(dp[2:4])  + float(0.5*45/3600)
latp = float(dp[:2]) * float(40/60) + float(0.5*30/3600)

lon = lonp + np.arange(80)*45/3600
lat = latp + np.arange(80)*30/3600

GeoTIFFの作成

GeoTIFFの作成はGDALを使いますが、（現時点では）Pythonicな記述ができないため、GDALがwrapされているrioxarray、fiona等を使ったほうがスマートです。今回はrioxarrayを使います。

GeoTIFFの圧縮は、 画像処理のど素人がまとめたpythonによるTIFF圧縮形式の比較 - 趣味で計算流砂水理 Computational Sediment Hydraulics for Fun Learning を参照頂きたいです。

import xarray as xr
import rioxarray

# y is North to South
ds = xr.Dataset({
                 'value': (['y','x'], x[::-1])                 
                }
                , coords={
                            'x': lon
                          , 'y': lat[::-1]
                         }
                    , attrs={'crs':'+init=epsg:' + str(4326)}
                )

dsp = ds['value']
d = dsp.rio.write_crs('epsg:4326', inplace=True)
out = dsp.rio.to_raster(f'geotiff\\{dp[:4]}.tif') #, compress='zstd')

GeoTIFFの結合

上記の手順で153個全てのshapefileをGeoTIFFに変換した後、取り扱いやすいように一つのファイルにまとめます。

方法は以下の2種類が考えられます。

VRTの作成

お馴染みのVRTを使う方法です。 個々のファイルが重い場合はこの方法が良いです。

from osgeo import gdal 

# gdal 2.1以降
my_vrt = gdal.BuildVRT('geotiff/output.vrt', glob.glob( 'geotiff/*.tif'), VRTNodata=65000, srcNodata=65000)
my_vrt = None

1つのGeoTIFFを作成

今回は個々のファイルが軽いので一つのファイルにします。 上記で作成したVRTから変換しています。

ds = xr.open_rasterio('geotiff/output.vrt')

d = ds.rio.write_crs('epsg:4326', inplace=True)
out = ds.rio.to_raster('geotiff\\AverageAnnualRainfall2010.tif', compress='zstd')

試してはいませんが、gdal_mergeでも出来るようです。 （参照：python - Stack geotiff images while retaining individual bands - Stack Overflow）

国土数値情報：平均年間雨量メッシュデータ

作成したGeoTIFFを図化すると以下のとおりです。

図化してみると非常に興味深いデータになりました。

  • 高降雨地帯は山脈沿いに集中している
  • 関東、関西の大都市圏は雨が少ない
  • 東北は太平洋側で雨が少ない
  • 屋久島は降雨量は5000mm級 ...etc

気象屋さんなら常識かもしれませんが。

気象庁｜報道発表資料のとおり、すでに「メッシュ平年値2020」が公開されています。 現時点では以下のpngのみですが、GML及びシェイプ形式が今後公開されるとのことです。 （気象庁｜メッシュ平年値図）

https://www.data.jma.go.jp/obd/stats/etrn/view/atlas/docs2020/png_large/prec/precipitation_13.png

大局的には2010年版と変わらないように思います。

GitHub

github.com

コードはこちら Jupyter Notebook Viewer

  • グラフの作成はHoloviewsを使用

pythonによる可視化はHoloviews一択 - 趣味で計算流砂水理 Computational Sediment Hydraulics for Fun Learning

  • はてブロへのhtmlグラフの埋め込みは以下を参照

はてブロにインタラクティブなグラフを埋め込む：外部htmlのとりこみ CDNとかiframeとか - 趣味で計算流砂水理 Computational Sediment Hydraulics for Fun Learning

備忘録：QGISのarea()と$areaの違い

R

たまに使うと忘れているので。

以下のサイトの一部をDeepLしただけです。

gis.stackexchange.com

楕円形と平面、どちらの面積も計算できます。

投影CRSで定義されたジオメトリの平面面積を知りたい場合は、area($geometry) 式を使用します。返される値は、その投影平面でジオメトリが囲む面積に相当します。

ジオメトリの楕円面積を知りたい場合は、$area 式を使用します。返り値は、楕円体 [2] 面の中でそのジオメトリ [1] が囲む面積に相当する。

どちらを計算するか、あるいはどちらを計算しなければならないかは、その値を知りたいかどうかによって決まります。

[1] 幾何学が投影型CRSで定義されている場合、幾何学の頂点は楕円に向かって逆投影され、結果として面積を計算する楕円形の幾何学は、それらの頂点の測地線の和で囲まれたものとなります。

[2] プロジェクトプロパティの計測セクションで定義された楕円体。

議事メモ：2022/11/5

R
  • 非構造格子の計算方法
    • よくわからないので勉強したい 。資料収集から。
  • 貯留関数のパラメータ同定
    • スパース推定を使えないか。
  • 浅水流モデルの高速化
    • ステンシル数が少なく並列化では限界がある。
    • dtを大きくとる、つまり、陰解法を用いるしかない。その場合、行列計算の高速化=並列化が鍵
  • 三次元地形データの活用方法
    • 高解像度の流れの計算が良い。高速化が鍵
  • オープンデータの話
    • メリットをマスにわかりやすく説明することが重要
  • AR,VRの話
    • 議論をマスに拡げるためには、想像により補間する部分を減らすために実世界と同じようにレンダリングすることは重要

画像処理のど素人がまとめたpythonによるTIFF圧縮形式の比較

R

画像処理なんて全くわからないのですが、必要になったのでまとめておきました。

TL;DR
  • TIFF画像作成時には多くの圧縮形式が準備されている
  • 今回の事例では、deflate、lzma、lzw、zstdによる圧縮でファイルサイズが30～50分の1になった。

はじめに

ジオメトリのメッシュデータをGeoTIFF化する場合は容量が非常に大きくなります。 ファイルサイズ縮小のため、GeoTIFFの圧縮は必要です。

圧縮方法は多くのアルゴリズムが提案されておりますので、いくつかの手法について概要と圧縮後のファイルサイズをまとめました。

TIFFの圧縮形式

pythonのラスタ処理用のライブラリrasterioでは 無圧縮を含めて、現時点で15種類の圧縮形式が準備されています。

参照：https://rasterio.readthedocs.io/en/latest/api/rasterio.enums.html#rasterio.enums.Compression

['none' ,'ccittfax4' ,'ccittfax5' ,'ccittrle' ,'deflate' 
        ,'jpeg' ,'jpeg2001' ,'lerc' ,'lerc_deflate' ,'lerc_zstd' 
        ,'lzma' ,'lzw' ,'packbits' ,'webp' ,'zstd']

これらの概要を簡単に下表にまとめました。

圧縮の仕組みまではよく理解していないです。

種類 概要
ccittfax4 FAXや同様の機器で電話回線を使用してモノクロ・イメージを転送するためのプロトコルとしてCCITT（国際電信電話諮問委員会）が開発
ccittfax5 同上
ccittrle modified Huffman run-length encoding : 同上
deflate zipと同じ形式
jpeg 代表的な画像形式。非可逆圧縮
jpeg2001 多分jpeg2000だと思う。jpegの高圧縮版で可逆圧縮も可能
lerc1,2 Limited Error Raster Compression (LERC)。 ラスターを複数のピクセル ブロックに分割する可逆または非可逆圧縮方式。
lerc_deflate lercとdeflateの組み合わせ？サポートしているソフトウェアが少ない。
lerc_zstd lercとzstdの組み合わせ？サポートしているソフトウェアが少ない。
lzma Lempel-Ziv-Markov chain-Algorithm. 7zと同じ形式
lzw Lempel–Ziv–Welch.圧縮効率と高速化の両面を追求している為、LZSSとハフマン符号化を組み合わせたDeflateアルゴリズム（LZHやZIP、PNGなどが採用）と比べると30%ほど圧縮効率が悪い
packbits 連長圧縮の改良版
webp 米Googleが開発しているオープンな静止画像フォーマット。非可逆圧縮でもアルファチャンネルを扱える。
zstd Facebookに所属しているYann Colletによって開発された可逆圧縮アルゴリズム。高速。

TIFFの圧縮サイズ

気象庁の1kmメッシュ解析雨量GPV：GRIB2形式をpython：hvplotでレンダリングする - 趣味で計算流砂水理 Computational Sediment Hydraulics for Fun Learningで使用した解析雨量のNetCDFファイルをGeoTIFFに変換して圧縮サイズの比較を行ってみます。

データは下図のとおりで、ピクセル数はlatitude: 3360longitude: 2560です。

以下のコードにより圧縮したGeoTIFFを作成しました。

import xarray as xr
import rioxarray

ds = xr.load_dataset('Z__C_RJTD_20180707100000_SRF_GPV_Ggis1km_Prr60lv_ANAL_grib2.nc')
ds = ds.rio.write_crs('EPSG:4326', inplace = True)

comp = ['none' ,'ccittfax4' ,'ccittfax5' ,'ccittrle' ,'deflate' 
        ,'jpeg' ,'jpeg2001' ,'lerc' ,'lerc_deflate' ,'lerc_zstd' 
        ,'lzma' ,'lzw' ,'packbits' ,'webp' ,'zstd']

for c in comp:
    out = ds.isel(time=0)['var0_1_200_surface'].rio.to_raster(c + '.tif', compress=c)

上記により出力されたGeoTIFFファイルの一覧は下表です。

  • deflate、lzma、lzw、zstdが高圧縮でファイルサイズが30～50分の1になっている。
  • その他はほとんど圧縮されていない。設定が必要か？
  • 今回のファイルはデータの大半がNaNと0のため、圧縮効率が良かったと思われる。
種類 ファイルサイズ 備考
none 32.832 MB
ccittfax4 0.026 MB ファイルエラー
ccittfax5 32.832 MB 無圧縮
ccittrle 0.026 MB ファイルエラー
deflate 0.418 MB
jpeg 0.002 MB ファイルエラー
jpeg2001 32.832 MB 無圧縮
lerc 32.832 MB 無圧縮
lerc_deflate 32.832 MB 無圧縮
lerc_zstd 32.832 MB 無圧縮
lzma 0.646 MB
lzw 1.602 MB
packbits 26.233 MB
webp 32.832 MB 無圧縮
zstd 0.372 MB

参考サイト

コーディングを意識したかなり丁寧なSIMPLE法：スタッガード格子版の離散化の解説

R

以前にSIMPLE法のプログラミングに関する以下の記事を書いたので離散化の記事を書いておきます。

SIMPLE法の記事は世の中に山ほどありますが、ソースコードとセットで書いている人は少ないのでコーディングをしたい方はぜひ読んで見て下さい。

TL;DR
  • SIMPLE法は離散化はかなり煩雑
  • SIMPLE法のポイントは運動方程式の線形化処理
  • 式展開の各段階で近似を行うため計算安定ため緩和係数が必須

はじめに

SIMPLE法は非圧縮性流体の計算方法として最も実務的な方法であると考えられますが、陰解法のため、陽解法と比較するとコーディングは難しいです。

私もそれほど頻繁に使用するわけではないので、久しぶりにプログラムを改良するとき等は結構気が重かったりします。 本記事では、自分自身の備忘録も兼ねて、プログラミングを意識して離散化についてまとめました。

なお、簡単のためx,yの二次元で記載しています。

運動方程式の離散化

次のx方向の運動方程式のみを記述しています。

\begin{align} &\phi = u \\ &\frac{\partial \phi}{\partial t}+\frac{\partial u \phi}{\partial x}+\frac{\partial v \phi}{\partial y} + \frac{1}{\rho} \frac{\partial P}{\partial x} - \frac{\partial }{\partial x} \left( \nu \dfrac{\partial \phi}{\partial x}\right) - \frac{\partial }{\partial y} \left( \nu \dfrac{\partial \phi}{\partial y}\right) = 0 \\ % \frac{\partial \phi}{\partial t}+\frac{\partial}{\partial x}\left(u \phi - \nu \dfrac{\partial \phi}{\partial x}\right) % +\frac{\partial }{\partial y}\left( v \phi - \nu \dfrac{\partial \phi}{\partial y}\right) % + \frac{1}{\rho} \frac{\partial P}{\partial x} % &= 0 \end{align}

スタッガード格子のため、コントロールボリュームは下図のように設けます。

まずは左辺第2,3項の移流項のみを抽出して離散化します。

\begin{align} \frac{\partial u \phi}{\partial x} + \frac{\partial v \phi}{\partial y} &= 0 \end{align}

として式展開を行います。

一次風上差分を適用するとx微分の項は次のように展開されます。

\begin{align} \frac{\partial u \phi}{\partial x} &= \frac{u^{n+1}_{Xp} \phi^{n+1}_{Xp} - u^{n+1}_{Xm} \phi^{n+1}_{Xm}}{\Delta x} \\ u^{n+1}_{Xp} &= \frac{u^{n+1}_{XP} + u^{n+1}_{C}}{2} \\ u^{n+1}_{Xm} &= \frac{u^{n+1}_{C} + u^{n+1}_{XM} }{2} \\ % \phi^{n+1}_{Xp} &= \left\{ \begin{array}{ll} \phi^{n+1}_{C} & (u^{n+1}_{Xp} \geq 0)\\ \phi^{n+1}_{XP} & (u^{n+1}_{Xp} \lt 0) \end{array} \right. \end{align}

プログラミング上は以下のように記述します。

\begin{align} u^{n+1}_{Xp} \phi^{n+1}_{Xp} = - max(-u^{n+1}_{Xp},0) \cdot \phi^{n+1}_{XP} + (max(-u^{n+1}_{Xp},0)+u^{n+1}_{Xp})\cdot \phi^{n+1}_{C} \end{align}

同様に、

\begin{align} u^{n+1}_{Xm} \phi^{n+1}_{Xm} = max(u^{n+1}_{Xm},0) \cdot \phi^{n+1}_{XM} - (max(u^{n+1}_{Xm},0)-u^{n+1}_{Xm})\cdot \phi^{n+1}_{C} \end{align}

y微分の項も同様に展開します。

\begin{align} \frac{\partial v \phi}{\partial y} &= \frac{v^{n+1}_{Yp} \phi^{n+1}_{Yp} - v^{n+1}_{Ym} \phi^{n+1}_{Ym}}{\Delta y} \\ v^{n+1}_{Yp} &= \frac{v^{n+1}_{YP} + v^{n+1}_{C}}{2} \\ v^{n+1}_{Ym} &= \frac{v^{n+1}_{C} + v^{n+1}_{YM} }{2} \\ \end{align}
\begin{align} v^{n+1}_{Yp} \phi^{n+1}_{Yp} &= - max(-v^{n+1}_{Yp},0) \cdot \phi^{n+1}_{YP} + (max(-v^{n+1}_{Yp},0)+v^{n+1}_{Yp})\cdot \phi^{n+1}_{C} \\ v^{n+1}_{Ym} \phi^{n+1}_{Ym} &= max(v^{n+1}_{Ym},0) \cdot \phi^{n+1}_{YM} - (max(v^{n+1}_{Ym},0)-v^{n+1}_{Ym})\cdot \phi^{n+1}_{C} \end{align}

これらを整理すると、

\begin{align} a_C \phi^{n+1}_{C} &= a_{XP} \phi^{n+1}_{XP} + a_{XM} \phi^{n+1}_{XM} + a_{YP} \phi^{n+1}_{YP} + a_{YM} \phi^{n+1}_{YM} \\ a_{XP}&= max(-u^{n+1}_{Xp},0) \Delta y \\ a_{XM}&= max(u^{n+1}_{Xm},0) \Delta y \\ a_{YP}&= max(-v^{n+1}_{Yp},0) \Delta x \\ a_{YM}&= max(v^{n+1}_{Ym},0) \Delta x \\ a_{C}&= a_{XP} + a_{XM} + a_{YP} + a_{YM} + (u^{n+1}_{Xp} - u^{n+1}_{Xm})\Delta y + (v^{n+1}_{Yp} - v^{n+1}_{Ym})\Delta x \end{align}

a_Cの最後の2項は連続式を用いると0となり、キャンセルされます。

次に移流項以外の項を含めて離散化します。離散式は次のとおりです。

\begin{align} &\frac{\partial \phi}{\partial t}+\frac{\partial u \phi}{\partial x}+\frac{\partial v \phi}{\partial y} + \frac{1}{\rho} \frac{\partial P}{\partial x} - \frac{\partial }{\partial x} \left( \nu \dfrac{\partial \phi}{\partial x}\right) - \frac{\partial }{\partial y} \left( \nu \dfrac{\partial \phi}{\partial y}\right) = 0 \\ &\frac{\phi^{n+1}_C - \phi^{n}_C}{\Delta t} +\frac{u^{n+1}_{Xp} \phi^{n+1}_{Xp} - u^{n+1}_{Xm} \phi^{n+1}_{Xm}}{\Delta x} +\frac{v^{n+1}_{Yp} \phi^{n+1}_{Yp} - v^{n+1}_{Ym} \phi^{n+1}_{Ym}}{\Delta y} \\ & + \frac{1}{\rho} \frac{P^{n+1}_{Xp} - P^{n+1}_{Xm}}{\Delta x} - \nu \frac{\phi^{n+1}_{XP} -2 \phi^{n+1}_{C} - \phi^{n+1}_{XM}}{\Delta x^2} - \nu \frac{\phi^{n+1}_{YP} -2 \phi^{n+1}_{C} - \phi^{n+1}_{YM}}{\Delta y^2} = 0 \end{align}

移流項と同様に式展開を行うと次のとおりとなります。

\begin{align} a_C \phi^{n+1}_{C} &= a_{XP} \phi^{n+1}_{XP} + a_{XM} \phi^{n+1}_{XM} + a_{YP} \phi^{n+1}_{YP} + a_{YM} \phi^{n+1}_{YM}+b - \frac{\Delta x \Delta y}{\rho} \frac{P^{n+1}_{Xp} - P^{n+1}_{Xm}}{\Delta x} \\ a_{XP}&= \left( max(-u^{n+1}_{Xp},0) - \frac{\nu}{\Delta x} \right)\Delta y \\ a_{XM}&= \left( max(u^{n+1}_{Xm},0) - \frac{\nu}{\Delta x} \right)\Delta y \\ a_{YP}&= \left( max(-v^{n+1}_{Yp},0) - \frac{\nu}{\Delta y} \right)\Delta x \\ a_{YM}&= \left( max(v^{n+1}_{Ym},0) - \frac{\nu}{\Delta y} \right)\Delta x \\ a_{C}&= a_{XP} + a_{XM} + a_{YP} + a_{YM} + \frac{\Delta x \Delta y}{\Delta t} + (u^{n+1}_{Xp} - u^{n+1}_{Xm})\Delta y + (v^{n+1}_{Yp} - v^{n+1}_{Ym})\Delta x \\ b &= \frac{\Delta x \Delta y}{\Delta t}\phi^{n}_C \end{align}

a_Cの最後の2項は連続式を用いると0となり、キャンセルされます。

SIMPLE法では定常解を取り扱うことが多いため、\Delta tを無限大にとって\Delta tを含む項を0とします。

連続式の離散化

同様に連続式も離散化します。

スタッガード格子のため、コントロールボリュームは下図のように設けます。添字の位置が運動方程式と異なる点に注意して下さい。

\begin{align} &\frac{\partial u }{\partial x}+\frac{\partial v }{\partial y} = 0 \\ &\frac{u^{n+1}_{Xp} - u^{n+1}_{Xm} }{\Delta x} +\frac{v^{n+1}_{Yp} - v^{n+1}_{Ym} }{\Delta y} \\ &(u^{n+1}_{Xp} - u^{n+1}_{Xm}) \Delta y +(v^{n+1}_{Yp} - v^{n+1}_{Ym}) \Delta x = 0 \end{align}

数値解法

手順1：仮の流速場の計算

前述の離散式の未知数は、u^{n+1},v^{n+1},P^{n+1}です。 これらは陽的に計算できないため、反復法によって求めます。

u,v,Pの仮の値をu^{*},v^{*},P^{*}とし、運動方程式に代入すると次のとおりになります。

\begin{align} \phi^{*} &= u^{*} \\ a_C \phi^{*}_{C} &= a_{XP} \phi^{*}_{XP} + a_{XM} \phi^{*}_{XM} + a_{YP} \phi^{*}_{YP} + a_{YM} \phi^{*}_{YM}+b - \frac{\Delta x \Delta y}{\rho} \frac{P^{*}_{Xp} - P^{*}_{Xm}}{\Delta x} \end{align}

手順1ではこの式よりu^{*},v^{*}について解きます。 この式はTDMAによる方法、反復計算による方法等で計算可能です。（前回のプログラムでは反復計算SOR法を使用）

この解法がSIMPLE法のポイントの一つになりますので詳細に説明します。

この計算では、仮の流速の反復計算の一つ前のステップの値\phi^{*old}から次のステップの値\phi^{*new}を求めています。 各項の係数aにはu^{*},v^{*}が含まれますが、上式の計算では反復計算の一つ前のステップu^{*old},v^{*old}のまま固定します。この線形化処理を行うことで上式が計算可能となります。

これらを用いて上式を書き直すと以下になります。

\begin{align} \phi^{*new} &= u^{*new} \\ a_C \phi^{*new}_{C} &= a_{XP} \phi^{*new}_{XP} + a_{XM} \phi^{*new}_{XM} + a_{YP} \phi^{*new}_{YP} + a_{YM} \phi^{*new}_{YM}+b - \frac{\Delta x \Delta y}{\rho} \frac{P^{*}_{Xp} - P^{*}_{Xm}}{\Delta x} \\ a_{XP}&= \left( max(-u^{*old}_{Xp},0) - \frac{\nu}{\Delta x} \right)\Delta y \\ a_{XM}&= \left( max(u^{*old}_{Xm},0) - \frac{\nu}{\Delta x} \right)\Delta y \\ a_{YP}&= \left( max(-v^{*old}_{Yp},0) - \frac{\nu}{\Delta y} \right)\Delta x \\ a_{YM}&= \left( max(v^{*old}_{Ym},0) - \frac{\nu}{\Delta y} \right)\Delta x \\ a_{C}&= a_{XP} + a_{XM} + a_{YP} + a_{YM} + \frac{\Delta x \Delta y}{\Delta t} + (u^{*old}_{Xp} - u^{*old}_{Xm})\Delta y + (v^{*old}_{Yp} - v^{*old}_{Ym})\Delta x \\ b &= \frac{\Delta x \Delta y}{\Delta t}\phi^{n}_C \end{align}

なお、係数a中のu^{*},v^{*}は反復計算で更新しますので、十分に収束していれば問題はないことになります。

手順2 : 圧力補正量の計算

次に、反復計算による補正量u^{'},v^{'},P^{'}を用いて、u^{n+1},v^{n+1},P^{n+1}を次のように定義します。

\begin{align} u^{n+1} &= u^{*}+ u^{'}\\ v^{n+1} &= v^{*}+ v^{'}\\ P^{n+1} &= P^{*}+ P^{'} \end{align}

これより、元の運動方程式と手順1の仮の値u^{*},v^{*},P^{*}を代入した運動方程式の差をとると次式が得られます。

\begin{align} \phi^{'} &= u^{'} \\ a_C \phi^{'}_{C} &= a_{XP} \phi^{'}_{XP} + a_{XM} \phi^{'}_{XM} + a_{YP} \phi^{'}_{YP} + a_{YM} \phi^{'}_{YM} - \frac{\Delta x \Delta y}{\rho} \frac{P^{'}_{Xp} - P^{'}_{Xm}}{\Delta x} \end{align}

補正量u^{'},v^{'}は十分に小さいと仮定して、それらを含む項を無視すると次式となります。

\begin{align} \phi^{'}_{C} &= - \frac{1}{a_C }\frac{\Delta x \Delta y}{\rho} \frac{P^{'}_{Xp} - P^{'}_{Xm}}{\Delta x} \end{align}

これが流速の補正量の式となります。

次にこの式を連続式に代入します。 まずは、わかりやすくするためにコントロールボリュームにあわせて添字を修正します。

\begin{align} u^{'}_{Xp} &= - \frac{1}{{a_C}|_{Xp}}\frac{\Delta x \Delta y}{\rho} \frac{P^{'}_{XP} - P^{'}_{C}}{\Delta x} \\ u^{'}_{Xm} &= - \frac{1}{{a_C}|_{Xm}}\frac{\Delta x \Delta y}{\rho} \frac{P^{'}_{C} - P^{'}_{XM}}{\Delta x} \\ v^{'}_{Yp} &= - \frac{1}{{a_C}|_{Yp}}\frac{\Delta x \Delta y}{\rho} \frac{P^{'}_{YP} - P^{'}_{C}}{\Delta y} \\ v^{'}_{Ym} &= - \frac{1}{{a_C}|_{Ym}}\frac{\Delta x \Delta y}{\rho} \frac{P^{'}_{C} - P^{'}_{YM}}{\Delta y} \end{align}

また、連続式は次のとおりに変形します。

\begin{align} (u^{n+1}_{Xp} - u^{n+1}_{Xm}) \Delta y +(v^{n+1}_{Yp} - v^{n+1}_{Ym}) \Delta x &= 0 \\ (u^{*}_{Xp} - u^{*}_{Xm})\Delta y + (v^{*}_{Yp} - v^{*}_{Ym})\Delta x + (u^{'}_{Xp} - u^{'}_{Xm})\Delta y + (v^{'}_{Yp} - v^{'}_{Ym})\Delta x &= 0 \end{align}

上式に前述の流速の補正量の式を代入すると次式が得られます。

\begin{align} A_C P^{'}_{C} &= A_{XP} P^{'}_{XP} + A_{XM} P^{'}_{XM} + A_{YP} P^{'}_{YP} + A_{YM} P^{'}_{YM}+B \\ A_{XP} &= \frac{1}{{a_C}|_{Xp}}\frac{\Delta y^2}{\rho} \\ A_{XM} &= \frac{1}{{a_C}|_{Xm}}\frac{\Delta y^2}{\rho} \\ A_{YP} &= \frac{1}{{a_C}|_{Yp}}\frac{\Delta x^2}{\rho} \\ A_{YM} &= \frac{1}{{a_C}|_{Ym}}\frac{\Delta x^2}{\rho} \\ A_{C}&= A_{XP} + A_{XM} + A_{YP} + A_{YM} \\ B &= - (u^{*}_{Xp} - u^{*}_{Xm})\Delta y - (v^{*}_{Yp} - v^{*}_{Ym})\Delta x \end{align}

手順2ではこの式よりP^{'}を解きます。 この式はTDMAによる方法、反復計算による方法等で計算可能です。（前回のプログラムでは反復計算SOR法を使用）

手順3 : 流速、圧力の更新

  • 手順2で求めたP^{'}により、 P^{*new} = P^{*old} + P^{'} として、P^{*}を更新する。
  • 手順2で求めたP^{'}から、手順2に示す流速の補正量の式よりu^{'},v^{'}を求める。
  • u^{'},v^{'}より、u^{*new} =u^{*old} + u^{'}、v^{*new} = v^{*old} + v^{'} として、u^{*},v^{*}を更新する。

手順4 : 反復

u^{'},v^{'},P^{'}が十分に小さくなるまで、手順1～3を繰り返す。

緩和係数の導入

SIMPLE法では各所で式の近似を行っているため、計算に不安定が生じやすいです。計算の安定性の観点から緩和係数を導入し、不足緩和を行います。

流速の不足緩和

全体：手順1～3の反復計算において不足緩和を行います。

具体的には反復計算の過程で得られた値\phi^{*new}を不足緩和係数\alpha_Uを用いて、次のように定義する。

\begin{align} &\phi^{*old} + \alpha_U(\phi^{*new} - \phi^{*old}) = \alpha_U \phi^{*new} + (1-\alpha_U) \phi^{*old} \rightarrow \phi^{*new} \end{align}

\phi^{*old}：反復計算の一つ前のステップでの値

よって、手順1の式は次のように変形できます。

\begin{align} \phi^{*new} &= u^{*new} \\ \dfrac{a_C}{\alpha_U} \phi^{*new}_{C} &= a_{XP} \phi^{*new}_{XP} + a_{XM} \phi^{*new}_{XM} + a_{YP} \phi^{*new}_{YP} + a_{YM} \phi^{*new}_{YM}+b - \frac{\Delta x \Delta y}{\rho} \frac{P^{*}_{Xp} - P^{*}_{Xm}}{\Delta x} \\ b &= \frac{\Delta x \Delta y}{\Delta t}\phi^{n}_C + (1-\alpha_U) \dfrac{a_C}{\alpha_U} \phi^{*old}_{C} \end{align}

圧力の不足緩和

流速と同様に不足緩和を行います。

手順3の圧力の更新時に不足緩和係数\alpha_Pを用いて、P^{*new} = P^{*old} + \alpha_P P^{'} とします。

補足：緩和係数の値

緩和係数は不足緩和のため1以下となる。対象とする流れ場によって最適値は異なるが、\alpha_U=0.5,\alpha_P=0.8程度と考えられている。

GitHub

GitHub - computational-sediment-hyd/SIMPLEalgorithm

nbviewer

参考書籍

バイブルです。読んだことが無い方はぜひ。

定番ですが。

参考サイト

時間ができたらコロケート版も書きます。

備忘録：python-openpyxlを使うときの注意-旧形式ファイル.xlsは使えません

R

pythonでexcelファイルを処理するときの最もポピュラーなライブラリopenpyxlを使用する場合、旧形式のexcelファイル（拡張子.xls）は使えないので、.xlsxまたは.xlsm形式に変換する必要があります。

大量のxlsファイルを変換する場合は、例えば以下のサイトの方法で、VBAを使って処理することが良いと思います。

【VBA】フォルダ内の.xlsを.xlsxに一括で変換する | 自恃ろぐ-jizilog.com-

備忘録：numpyでインデックス用の連番を振る

R

数値計算ではインプットファイルなどでインデックス用の連番を振ることよくあります。

テクニックのメモです。

配列を繰り返す

配列0～4を3回繰り返したインデックスの作成方法です。

import numpy as np

x = range(5)
print(list(x))
# [0, 1, 2, 3, 4]

i = np.tile(x, 3)
print(i)
# [0 1 2 3 4 0 1 2 3 4 0 1 2 3 4]

配列の要素を繰り返す

配列0～4の要素を各3回繰り返したインデックスの作成方法です。

i = np.repeat(x, 3)   
print(i)
# [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4]