地下水の熱輸送計算 その2

1D の続きです。

 地下水の熱輸送算結果は以下の通り。50度の地下水を間隙率10％の地盤に1m/sec で注入しました。図の見方は先日と同じです。

 

 最初、100秒程度では、境界付近しか温まらない計算結果が出ました。遅すぎて驚きましたが、そういうものかもしれないと気を取り直し、100,000秒計算。でも遅い。

dtを1000秒にしていたのであまりにも飛ばしすぎたかと細かくしましたが、結果は同じ（同じなのも意外ですが）。うーん、あやしい。

段ボールに戻り、再度本探し。次に見つけたのがコチラ↓。

David E. Radcliffe「Soil Physics with HYDRUS: Modeling and Applications」

式4.56が移流拡散方程式です（ちなみに有機物量や不飽和にも対応しています）。
よく見ると、利用していた登坂「地圏水循環の数理」式8.35に比べ輸送項の係数が異なります。
どちらが正しいのか？と思いながら単位を確認すると、「Soil Physics with HYDRUS」があっていました。内部エネルギーのT-T0の部分が時間項の差分化部分と捉えると、係数は残る比熱に密度をかけて体積熱容量になっているので、全体でJ/m3/sになる必要があります。そうすると、流速に比熱を乗ずるのは納得。誤植でしょうか？というか、式見て気づけよ、私！

二重に誤っていたのが流速。実流速かと思いきや、ダルシー流速と書かれています。これはどちらの図書にも導出が記載されていませんので、あとでチェックですね。

気を取り直して再計算。

流速1m/s、間隙率100％。濃度と同じ、流体のみの計算です。境界条件も濃度と同じように時系列変化させました。

OKです。定常で均一に近くなるというのも正しかったようです。次は間隙率10％に。

熱輸送速度をみると単純に10倍にはなっていませんが、速くはなりますね。

ダルシー流速1/10とみなして0.1m/s、間隙率10％。

 

遅い気がしますね。実流速、ダルシー流速のチェックをしておかないと心底納得できないかも。 

次は2D断面。
コチラ↓の方と同様ですが、左端の流体部分は地盤を考慮して係数を修正しています。
http://penguinitis.g1.xrea.com/study/FiPy/FiPy.html

左端に断層をイメージした高透水帯を設け、下から温泉を注入。地表部（上端）の境界はこれまで同様に適当な温度勾配を設定しただけです。その他も適当なので得られた値に厳密性を期待していませんが、なんとなくイメージが湧いてきました。

 今回はココまで。地下水の熱輸送計算にFiPyを利用てみました。個人的にはHydrus 1D よりも拡張性が高いと感じています。PHREEQC の COM 版とあわせると面白そうですが、それはPHREEQC 単体の方が簡単ですね。1～2次元の計算で候補といったところでしょうか。簡単な試算か、土石流や河川、浸透流の結果を取り込んで計算する使い方が主体になりそうです。3次元なら Dtransu を修正した方が良いでしょう。

わからないところや設定方法を詰めないといけないところ（境界条件）も出てきました。あらためて調べておき、納得できる答えを出せるようにしましょう。

 

地下水の熱輸送計算

この週末で、寝かしていた地下水の熱輸送計算を片付けておこうと思い立ちました。

利用するのは FiPy。
1次元だと EXCEL でもよかったのですが、いまさら感があります（余談ですが、ここ数年で取り扱うデータ数が飛躍的に増加し、 EXCEL はベンチウォーマーになっています）。

伝熱計算では、以前見かけたコチラ↓のサイトに簡単なコードが紹介されています。 http://penguinitis.g1.xrea.com/study/FiPy/FiPy.html

１次元でも2Dメッシュを使った方がビジュアル的に理解しやすいので、そのまま2Dメッシュを使わせていただくことにしました。

地表流などでの熱の移流拡散、定常・非定常はそのままでOKですが、地下水など媒質中を流れる流体の熱の移流拡散は、各項の係数に修正が必要となります。その修正は簡単なのですが、あらためて調べてみるとあまり紹介されていないですね。今回は、段ボールの中から最初に見つけた以下の図書を参考にしました。

登坂「地圏水循環の数理」

で、動かしたのですが、定常の結果がほぼ一定値。左側（流入側）の境界条件に入力した値より1K低い値が連続します。
左側の境界条件は温度固定なので一定値でもOKです。が、右側（流出側）の境界条件は温度固定ではなく勾配を適当に入れていましたので、ある範囲で温度勾配がつくように予想していました。これでは結果があっているのか計算ミスをしたのかわからず、迷ってしまいます。

ひとまず係数をもとに戻し、理解しやすい表流水内の濃度の1D移流拡散でスクリプトを整理することにしました。結果は以下の通りです。
上図は最終ステップの濃度分布。2Dメッシュなので色付きで表示してくれます。下左図は縦軸にTimeStepをとって時系列変化を表現したモノ（STIVと同じ絵です。https://phreeqc.blogspot.com/2019/12/river.html）。下の参照リンクと異なり転置していません。
あと、time step も1秒と飛ばしているので、得られた濃度に厳密な意味は期待しません。が、与えた流速と整合する絵が得られましたので計算は正常に回ったと判断します。

ここで悩んだのは左端の境界条件。時系列で変化する濃度（10秒まで0, 50秒まで50, 60秒まで, 100, 以降50）を与えたかったのですが、その方法がわかりませんでした。が、検索すると同じことを聞かれている方がいらっしゃいました。コード付きの回答があります。ありがたい。
https://stackoverflow.com/questions/56237476/non-steady-diffusion-advection-equation-with-time-dependent-dirichlet-boundary-c

ただ、このコードでは境界条件が更新されなかったので、左端の境界条件設定文を for loop 内に入れたら期待通りに動きました。

スクリプトは組めたので、再度、熱輸送計算へ。
続きは後日。

 

 

Newmark + Grid Search

先日から触っている簡易なニューマーク法は、Newmark's sliding-block model とか、Newmark (rigid-block) analysis などと呼ぶそうです。

その変位は Newmark displacement で統一されているようですね。国内ではどのように呼ぶのでしょうか？

この簡易な手法、 Pythonなどの言語（+GIS） と相性が良いと思います。
今回、c・φ、水位などの取り得る範囲はわかるけれども妥当性がつかめないパラメータは、すべてグリッドサーチにかけました。今までならパラメータ設定や根拠付けに時間をかけていたところでしょう（そしてイマイチな精度を出していたことでしょう）。機械学習を経験してからは、合えば良いという割り切りを所々で使うようになりました。

今後、こういった考え方は Python や機械学習を利用する若い世代が育つにつれ拡大するかもしれません。そうすると、今まで適用性が狭く精度の低かった分析が容易に広く使われだすかもしれません。発想の転換が可能性を広げます。

人も技術も順調に育ってほしいですね。期待しています。

地震時崩壊

10年ぶりの余震、幸い被害は以前に比べ小さいようでした。

ちょうどプロと話をしながら地震時崩壊について試算していましたが、ここまで崩壊が少ないと何をやっているのだろうと思いますね。計算する前からまったく合わないことがわかってしまいます。何か見逃しているのでしょう。

地震時の地すべりや崩壊を扱う研究は多くあります。土研さんや国総研さんは以下の資料を出されています。両方とも手を動かした結果、国総研さんの手法は精度が良くないもののある程度正しそう、というのが個人的印象です。

既存地すべり地形における地震時地すべり発生危険度評価手法に関する研究（土研資料）
https://www.pwri.go.jp/team/niigata/dokensiryo4204_web.pdf
https://phreeqc.blogspot.com/2012/11/blog-post_18.html

地震による斜面崩壊危険度評価手法に関する研究（国総研資料）
http://www.nilim.go.jp/lab/bcg/siryou/tnn/tnn0204.htm

海外ではコチラ。今月の文献です。
Near-real-time prompt assessment for regional earthquake-induced landslides using recorded ground motions
https://www.sciencedirect.com/science/article/pii/S0098300421000236
Newmark法 が使われています。プロによれば波形補間部分の技術論はイマイチ（結果もイマイチ）なのですが、チャレンジ精神には拍手です。
既存文献を追いましたが、即時予測を目的とした簡易な Newmark 法は古くから研究されているようですね。例えばUSGSさん。
https://pubs.er.usgs.gov/publication/70044002
これ、１地震だけ手を動かしてみましたが、国総研さんの手法よりは良い結果が出ました。面白い。

いずれも精度は高いといえないものの、簡易で解釈性も高い手法ばかりです。機械学習で精度を追求する一方、こういった古くからある手法での試算を改良し、まだ見逃している項目を探す必要があるのでしょう。

matplotlib

matplotを利用する際、毎回どこかの例をコピペするので、ココにまとめておきます。
今後、追加していきましょう。

import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
from matplotlib.font_manager import FontProperties
fp=FontProperties(fname='c:/WINDOWS/Fonts/msgothic.ttc',size = 10)
#fp=FontProperties(fname='./font/TakaoGothic.ttf')
#import japanize_matplotlib
#import seaborn as sns
#sns.set(font="IPAexGothic") 

#####################################################
#ヒストグラム
#CDF（2軸）
#凡例は外側に配置、透明度0
#####################################################

def plot_hist(col,SaveFName):
    
    data=df[col]
    fig, ax1 = plt.subplots()
    ax2 = ax1.twinx()

    # histogram
    n, bins, patches = ax1.hist(data,
                                bins=200,
                                #density=True,
                                #cumulative=True,
                                range=[-100,100],
                                alpha=0.7,
                                label='頻度')
    # another y-axis
    y2 = np.add.accumulate(n) / n.sum()
    x2 = np.convolve(bins, np.ones(2) / 2, mode="same")[1:]

    lines = ax2.plot(x2, y2, ls='-',
                     color='r',
                     #marker='o',
                     label='累積分布')
    # set axis
    ax1.set_xlim(-100, 100)
    ax1.set_ylim(0, ylim)
    ax2.set_ylim(0, 1.05)
    #ax2.yaxis.set_major_locator(ticker.MultipleLocator(0.2))
    ax1.grid(visible=True)
    ax1.set_axisbelow(True)

    # legend
    handler1, label1 = ax1.get_legend_handles_labels()
    handler2, label2 = ax2.get_legend_handles_labels()
    ax1.legend(handler1 + handler2,
               label1 + label2,
               loc=2,
               shadow=True,
               #borderaxespad=0.
               prop=fp
              )

    ax1.set_xlabel(col, fontproperties=fp)
    ax1.set_ylabel('頻度', fontproperties=fp)
    ax2.set_ylabel('累積確率', fontproperties=fp)
    plt.savefig('./output/'+SaveFName+'.png',
                format='png',
                dpi=300)
    plt.show()


#####################################################
#棒グラフ
#CDF（2軸）
#凡例は外側に配置、透明度0
#####################################################

fig, ax1 = plt.subplots(figsize=(6,3)) 
ax2 = ax1.twinx()
i=0
        
for col, data in df.iteritems():
    data_sum=data.sum()
    data1=data.div(data_sum)
    data_cum = data1.cumsum()
    
    bar= ax1.bar(df.index+0.4*i,
                         data,
                         width=0.4,
                         align='edge',
                         label=str(col))
    line = ax2.plot(df.index+0.5,
                data_cum,
                ls='--',
                #marker='o',
                #color='r'
                label='累積'+str(col))
    i=i+1
                         
# set axis
ax1.set_xlim(10, 100)
ax1.set_ylim(0, 105)
ax2.set_ylim(0, 1.05)
ax1.grid(visible=True)
ax1.set_axisbelow(True)

# legend
handler1, label1 = ax1.get_legend_handles_labels()
handler2, label2 = ax2.get_legend_handles_labels()
ax1.legend(handler1 + handler2,
                       label1 + label2,
                       bbox_to_anchor=(1.15, 1.02),
                       loc='upper left',
                       shadow=True,
                       #borderaxespad=0.
                       prop=fp
                      ).get_frame().set_alpha(1.0)

ax1.set_xlabel('aaa', fontproperties=fp)
ax1.set_ylabel('bbb', fontproperties=fp)
ax2.set_ylabel('ccc', fontproperties=fp)

plt.savefig('./output/ddd.png',
            format='png',
            bbox_inches='tight',
            dpi=300)


#####################################################
#折れ線（マーカーのみの散布図）
#移動平均
#dfのカラムを順にプロット
#散布図と移動平均の色を合わせる
#####################################################

def plot_rolling(df, SaveFName):
    fig, ax1 = plt.subplots()

    for i in range(df.shape[1]):
        lines = ax1.plot(df.index,df.iloc[:,i],
                         ls=' ',
                         color='C'+str(i),
                         marker='o',
                         label=df.columns[i])

    for i in range(df_rate.shape[1]):
        s=df.iloc[:,i].rolling(window=10).mean()
        lines = ax1.plot(df.index-4.5,s, 
                         ls='-',
                         color='C'+str(i),
                         #marker='o',
                         label='移動平均'+str(df.columns[i]))

    # set axis
    ax1.set_xlim(0, 50)
    ax1.set_ylim(0, 0.05)
    ax1.grid(visible=True)
    ax1.set_axisbelow(True)

    # legend
    ax1.legend(bbox_to_anchor=(1.05, 1), loc='upper left',prop=fp)

    ax1.set_xlabel('aaa', fontproperties=fp)
    ax1.set_ylabel('bbb', fontproperties=fp)

    plt.savefig('./output/'+SaveFName+'.png',
                    format='png',
                    bbox_inches='tight',
                    dpi=300)