継続時間、再現期間と崩壊

2011年、台風12号による奈良県の深層崩壊では、”雨”についていくつかの指摘がありました。

その中で、2軸指標（長期雨量、短期雨量、土壌雨量指数などの組み合わせ）を使うという動きがあります。「技術者に必要な斜面崩壊の知識」にも簡単ですが掲載されていました。p215のグラフです。

それをもっとわかりやすくしたのが、次ページの降雨継続時間と再現期間での議論。これは面白いですね。

今年、初めて自分で手を動かして確率降水量を計算した際に、継続時間の区切りなどで引っかかったことがあります。この本では、それを意識することなく、継続時間毎に再現期間を計算しプロットしてありました。これにより、一連の雨の中で、どの継続時間がもっとも斜面にとって不慣れな雨だったかを視覚的にわかるように工夫されています。砂防分野ではスタンダードな手法なのでしょうか？

3章では確率降水量の基礎、再現期間と表層崩壊の話がまとめられています。このように、継続時間毎の再現機関を比較できると、どの継続時間の雨が斜面崩壊に起因したかを推定しやすくなります。当然、概略の域を超えませんが、傾向をつかむには良い手法だと思います。また、これらが土研さんが提供しているEXCELシートで計算できることも紹介されています。
http://www.pwri.go.jp/jpn/seika/amedas/top.htm

実際に触ってみましたが、非常に手軽ですね。アメダスに限定されますが、データを集めて期間毎に整理する必要がありません。国土技術センターさんの水文統計ユーティティとあわせて持っておけば、多くの場合、事足りるでしょう。

台風12号の深層崩壊では、72時間の再現期間が約200年のようです。その前に、48時間が150年ですので、雨量だけではどのあたりで崩壊したかが分かりません（振動センサーなど、他の観測結果で分かっていますが）。しかし、傾向として、短期より長期の雨量が崩壊に対し危険ということはすぐにわかります。すばらしい。

一方、誘因としての再現期間200年では、斜面が安定していたと考えられる数千年（数千年の根拠は書かれていません）で壊れていないのが説明できないとも。従って、クリープなどの素因の変化があるはずと展開されています。ま、この辺は地質屋さんが根拠を持ってこないといけませんね。

古いデータを使われている箇所もありますので、この手法自体、昔からあるのでしょう（私が知らなかっただけです）。もう少し、読み進めましょう。

斜面の動水勾配

大掃除も終わり、朝から家でゴロゴロ。
振り返ってみると、これだけゆっくりしたのは、この一年でも久々です。年末ですね。

冬休み用に2冊本を買っていました。そのうち、1冊目を読みはじめました。

飯田智之「技術者に必要な斜面崩壊の知識」

引っかかったところが1か所。p142の動水勾配について。
いえ、私の基礎知識がなかったわけなのですが、最初はそれに気が付きませんでした。

非常に基礎的かつ単純な話です。
斜面で基盤岩に沿った流れがある場合、その動水勾配をsinθで表現されています。理由は分かりますが、tanθでは間違いなの？と思いつつ調べてみますと、やはりsinθでないといけないようです。こういった問題は、学生のころに習うようですね。以下は新潟大学の講義資料です。
http://geotech.eng.niigata-u.ac.jp/lecture/enshu12/ans121107.pdf

ΔLを水平距離と考え，動水勾配を tan 20°とした間違いが多く，正答は7名のみでした。
この問題では，流線がすべて斜面と平行になっています。動水勾配を求める際の透水距離ΔLは流線に沿った長さなので，ここでは水平距離ではなく斜面長になります。
したがって，動水勾配は浸潤面（斜面）勾配のtanではなくsinの値になります。
なお，動水勾配が正しいが，透水断面積を 2.8 m2としていた解答も目に付きました。

見事に動水勾配を間違え、断面積にcosをかけるのを間違えました。情けないですが、実力です。独学ですと、こういった基礎力がぽっかり抜けていることがあるんですよね。自覚はしているのですが、なかなか見つけられません。今回は良い機会でした。

もう少し読み進めましょう。

少し早い仕事納め

今年の仕事納めは例年より早く、26日でした。
その後はメールでのやり取りが続き、メールの来なくなった今日が本当の仕事納め。現場で締めなかったのでやや拍子抜けしていますが、ま、こういう”危険”のない年末もたまには良いでしょう。

今年もいろいろ経験させてもらいました。
前半は現場とシミュレーションが半々、夏から秋にかけて現場が続き、最後にシミュレーションといった、変則的な年でした。

中期的な目標「多くの現場で経験をつむ」に対しては、そこそこの年でした。もっと多くの現場を経験しないとダメですが、なかなか携わることのない分野、スケールの現場にかかわることができました。個人的に良い経験だったと思います。
また、意図せずシミュレーションの種類や幅も増えました。移流分散や土砂移動のシミュレーションコードを改良し、満足いくレベルで実施可能になったのは今年の収穫でしょう。

短期目標としていた技術士試験はダメでしたね。来年も同じ目標に設定します。ただ、選択問題減に対して対処する方法はわかるのですが、それに対するモチベーションが上がっていません。最初は無理に手を動かして様子を見てみましょう。


振り返ってみると、人との縁で、思わぬチームに混ぜていただいたり、声をかけていただいたり。逆にギクシャクしたり。営業の方のように、うまく対応できるスキルがあれば良いのですが、嘆いていても仕方ありません。めげずに来年も頑張りましょう。

libfcoremdd.dll が見つからなかったため、このアプリケーションを開始できませんでした。

後輩から、以下の報告がありました。

「”アプリケーションを正しく起動できませんでした(0xc000007b)。”というエラーが出ます。」

私がコンパイルした exe で計算しようとし、このエラーが出た模様です。
似たようなことがあったなぁと思いつつ、このブログを検索すると、過去に同じエラーを解決していました。
http://phreeqc.blogspot.jp/2012/08/0xc000007b.html

しかし、今回は同じように Redistributable Libraries x64 をいれてもらっても、解決しませんでした。私の使用しているPCより新しい XEON が入っているので、そのせいか？と思いつつ、別の PC で試してもらうと、今度は別のエラー。

”libfcoremdd.dll が見つからなかったため、このアプリケーションを開始できませんでした。”


意味が分からずいろいろ調べ、やっと原因特定。
32bitでコンパイルしていました。それが x64 の dll を見に行っていたため、0xc000007b のエラーが出たわけです。
しかもデバッグ版をわたしていました。libfcoremdd.dll が出たところですぐに気づくべきでしたね。まだまだです。

すまない、後輩君。来年は、もっと頑張ります。

技術者の階段

御飯を食べながら録画していたTVを聞いていると、おおぉ、というセリフが。

リーガルハイ2で伊東四朗さん（被告側）が言われたセリフ。

そもそも才能なんてものはな、自分で掘り起こして見つけるものなんだ。
俺だって天才なんかじゃない。
誰よりも必死で働き、階段を１つ１つ、踏みしめてきただけだ。
振り向いたら、誰もついてきてない。
怠けた連中がふもとでこうつぶやく。
「あいつは天才だから」
冗談じゃない！

正論です。
が、サラリーマン技術者は言えないセリフです。基準はあくまで努力している技術者と、していない方の平均、あるいは上司の技術レベルや価値観です。

先日の忘年会で、後輩の大残業、ストレス太りの話が出ました。「上司が仕事をしないので、仕方なく上司の仕事をやっている」と嘆いていたのを常に聞いていましたし、仲間内ではその上司が仕事をしないことで有名でした。実際、昨年度は私の1/10以下の生産でしたので、十分予想できたことです。しかも、太りだしたことを言い出したのが仕事を後輩に頼んでいるその上司自信。後輩は何も言えず。
で、思わず、その上司に「もう少し仕事がんばってもらえたら、」と言ったら、激高されました。コンプレックス ＆ NG ワードだったのでしょうか。プライドを大きく傷つけてしまったようです。その後はコンコンと、自分は頑張っていることを主張されていました。謝ったのですが、それでも怒り心頭のようでした。別の上司や後輩はその上司より仕事を多く抱えていたのですが、もうそれ以上話すのはやめました。そういえば、昨年度も別の方に突っ込まれて、憤慨されていたようです。その位置における、それなりの言い分はあるのでしょう。

将来、他の方が階段のどの位置にいて、自分がどこにいるということを正しく理解するとともに、そこで生まれた差が何だったのかを判断し、努力し続けることが、技術者として必要だなあと感じたセリフと出来事でした。


もう一つ。原告側のセリフ。

すぐに追い抜いてやる、王様の椅子は俺がもらう。
あんたのより遙かにどデカいピラミッド作ってやるよ！

コチラはセリフ自体に別に共感したわけではありません。ピラミッドというキーワードに引っかかりました。解釈は異なるかもしれませんが、はるかにでかいピラミッドを作ろうとすると、基礎の幅を大きくする必要があるなあと、いまさらですが感じました。
地質屋さんは地質だけをやっていてプロになれるわけではありません。数学・物理・化学の基礎力、工学など他分野の知識や表現技術、あとは現場での経験を大事に、幅広く身に着ける必要があります。基礎を固め広げ続けることと、コア技術を伸ばすこと、両者が技術者には必要なのでしょう。

良い気分転換になりました。

========================================================
2013/12/29追記
会社に行くと、激怒した上司とは別の上司が出社されていました。明日も現場へ行くとのこと。私もまだまだ努力が足りません。

2013/12/31追記
ここまではひどくないですが、訳アリ感含めて似ていますね。ヨイショしとけばいいんでしょうけど。
http://detail.chiebukuro.yahoo.co.jp/qa/question_detail/q12106710507

Autodesk Geotechnical Module

地質を3次元で可視化する場合、その根拠となるボーリングも同時に表示します。

表現手法は様々ですが、個人的には MVS のように、シリンダーに球がつく表現法が好きです。球は濃度、N値、透水係数など目的によって変えますが、基本的にはスカラーの表現です。

GEORAMA のような簡易柱状図表示は、見る方向が限定されますし、細かいN値などが読めなくなるので、あまり好きではありません。
そのため、今回、MVSを使わずに infraworks へデータを持っていくには、どのような形にすれば見栄えがよいか？と考えていました。

思い付いたのが Geotechnical Module。以前、GEORAMAの洗礼が嫌になり、一度触ってみたことがあるのですが、使用法がよくわからなかったので放り投げていました。
今回、ヘルプを読んでみますと、そのコンセプトは GEORAMA と同じでした。ボーリングから地層のサーフェス作成、およびそれらのデータ管理といった内容です。操作法は動画でも詳しく説明されています。
http://kbs.keynetix.com/category/term/488/488

無償ですので、機能は GEORAMA に勝てません。ボーリングの地層境界を直線で結んでTINを作るといった簡易的な面の作成はできますが、複雑な形状は手間がかかりそうです。が、柱状図は地質による色分けシリンダーになっています。これがSolidであれば問題なく infraworks に移行できそうです。

これを確認しようと、csv を読み込むまで実施しましたが、そこでやめました。データサーバーとライセンスマネージャーを別途インストールする必要があったのです（最後まで ヘルプを読んでからやればよかったのですが）。

とりあえず、今日はここまで。別の方法がなければ続きを行いましょう。

==========================================================
2013.12.28追記
結局、MVS でボーリングを表示してしまいました。
来年、Infraworks を触る場合には、MVS から CAD データに書き出しましょう。

地質の可視化

先週より、シミュレーションコードの改良と、3次元可視化の作業を行っています。

コードは簡潔で分かりやすい。計算が軽いのも手伝い、私のような素人でも非常に読みやすく、簡単に改良できました。こういったコードをかけるようになりたいですね。

可視化のほうはイマイチ。先のシミュレーションで使用した3次元地質モデルを甘めに、さらに広域で作成しています。が、面倒。しかも one off （これが気に入らない）。
2次利用とはいえ、可視化だけのためにモデルを作成するのは初めてです。今までは、計算という目的や、正しい2次元モデルを、3次元を利用して作るといった目的がありました。今回は見せるためとはいえ、ただ作るためだけに作るといった感があります。興味はそれほどありません。
CG専門の方が作成されたら良いと思うのですが、調査結果のない箇所も含めた地下のモデル化というのは、やはり地質屋でないと難しいようです。逆に言えば、地質屋がこういった見せる努力を行ってこなかったので、一般の方も見える部分だけに意識が向き、見えない地下に意識が向かなかないのかもしれません。

広域のモデル化は重たいですね。空中写真を読み込むだけでも時間がかかります。
本来、こういった作業はCADではなく、可視化に特化したツールを使わないといけないのでしょう。地形の可視化＋地質表示だけなら、ツールは多くあります（私が使えるのは限られていますが）。
エンドユーザーレベルでもいくつかありますよね。例えば、Google Earth。精度は粗いのですが、携帯端末でも高速で動きます。ただ、写真は被災後、DEMは被災前なんてことがほとんどですので、被災地などの可視化といった目的には不向きです。ま、それが目的で作られているわけではないので贅沢は言えません。

ArcGIS で可視化が可能だろうと思い、モデル（dwg）を読んでみましたが、ソリッドは表示できませんでした。shp か何か、読める形にする必要があるのでしょう。
ん？取り込んでしまえば、 MODFLOW に持って行って、結果を MVS で可視化するなんてことも可能では？試してみる価値はありそうですね。

GEORAMA ユーザー

構造物を3次元ソリッドで作成していただこうと、他分野の方と話をしていました。

席に来られるなりPCの画面をのぞかれ、「Civil3D ですか」と。
グラボ、CPU の名称を話しても通じました。ケースも珍しそうに見られていました。GEORAMA も御存知でした。話を聞くと、Civil3D や InfraWorks をよく使用されているそうで、ユーザー会にも積極的に参加されているようでした。

その方のお話ですと、GEORAMAユーザーはレアとのこと。
ま、そうでしょう。BIMやCIMなど、設計者による3次元化は始まっています。が、上物が主体です。地質屋さんががんばって3次元化を進めているといった話も聞きませんし、必要だという設計者の声も聞きません。

ただ、GEORAMAも、それしかないので使っているだけですし、おすすめできるソフトではありません（先日もエラーで立ち上がりすらしませんでした）。
ですが、3次元で地質分布を組み立てていく方が、間違いが少なくて済みます。断面毎の整合性や分布のなめらかさをPCにチェックさせた方が楽ですし、自分のミスに気づきやすくなります。とにかく、手を動かさないと気付かないことだと思います。

今週は、大手さんが調査された結果を使って、シミュレーション用の3次元地質モデルを作成していました。ですが、3次元で地質を組み立てている = 他社の地質屋さんのミスを手直ししている、といったような状況でした。横断で断層のずれを考慮しているが縦断で断層が表示されていない and 断面のズレが考慮されていない、断面交差位置で崩土の厚さが違う、平面では1つのすべりブロックであるが、断面では途中から段差が生じるなど。それ以外にも基本的なこととして、ボーリングの位置が読み込んだXMLと平面で違う、平面に座標を持たせていない＋トンボがない=大体でしかプロットできない、等々。気付かれていないのでしょうね。


Civil3D2013 ＋ GEORAMA2013 （いまだにGEORAMAの2014が出ていません。どうなっているのでしょう、CTCさん）でモデルを作った後の可視化については MVS を考えていたのですが、その方曰く、「InfraWorks のほうがきれいかも」。
データを移す手間はほとんどないそうですので、手を付けてみますか。

深層崩壊と隆起量

先日、「深層崩壊」を見てきました。

といっても、崩壊の発生は数年前であり、対策は進んでいます。
今回は崩壊面の位置と状況を把握することがメインでした。LPだけでも読めるのですが、やはり、現場を見ると情報量・解釈のしやすさが違います。

７か所ほど回ったのですが、どこも一様ではありませんでした。流れ盤である箇所、受け盤である箇所、１回で落ちている箇所、２回以上の箇所、砂岩主体の箇所、頁岩が多い箇所、岩盤崩壊がメインの箇所、表層の土砂崩壊も多く含まれている箇所、天然ダムの有無など様々です。

ただ、共通しているのが地形。
山の形状は隆起量の大きい地域と似ていましたので、ここもそうなんだろうなあと想像しながら歩いていました。
帰ってから調べてみると、やはりそのようでした。隆起量の大きな地域は、まさに今、山が高くなろうとしていると同時に、不安定になった箇所は削られ険しくなっている地域と言えます。
国交省さんの発表資料でも、深層崩壊と隆起量の関係をまとめられていました。
平成24年度 国土技術政策総合研究所講演会「深層崩壊～その実態と対応～」
http://www.nilim.go.jp/lab/bbg/kouenkai/kouenkai2012/pdf/121204_6.pdf

「深層崩壊」という用語はよく使用されていますが、いまだにその定義はあいまいなようです。崩壊土砂量で定義していることもあるようです。砂防学会によれば「山地および丘陵地の斜面の一部が表土層（風化の進んだ層）のみならずその下の基盤を含んで崩壊する現象」だそうです。
http://www.jsece.or.jp/news/2012/teigen20120402.pdf
いずれも、規模の大きな崩壊や地すべりを指すことが多いようです。


移動中に、車道からとった写真です。三角形のところに,山中の天然ダムが見えます。大きいですよね。

Phast４windows

PHAST の最新版を DL しようと、USGS の HP をのぞいてみました。
http://wwwbrr.cr.usgs.gov/projects/GWC_coupled/phast/

実際に使用するのは久しぶりでしたが、気になるので時々のぞいていました。
64bit 版あり、GUIありと、定期的に更新されています。

で、以前から気になっていた「 Phast４windows 」というプレを DL。
解凍してみると、なんと WPhast の最新版でした。WPhastは参考書についていたプレですが、使い勝手がよく、気に入っていました。ですが、最新版を DL する箇所が分からず、ずっと 2008年の Ver.0.6.1 を使用していました。

使ってみると、WPhast とほぼ同じインターフェース。非常に簡単です。浸透流と PHREEQC を御存知の方であれば、すぐに使えると思います。

Arc の shp も読めるんですね。shpが読めるなら、地層のモデル化もできそうですね。

唯一、 parallel の計算だけは動きませんでした。parallel は浸透流と反応を平行で計算するという意味でしょうか？多コアでしょうか？これはマニュを読まないとダメでしょうね

ポストは Model_Viewer。これも変わっていません。新規作成し、h5を取り込むだけ。
http://phreeqc.blogspot.jp/2010/10/as.html＜こんな感じです。



で、肝心のPRBsの結果はというと、イマイチ。
うーん、仕方ないですね。

透過反応壁とPHREEQC その2

結果が微妙に異なるのは、使用していた熱力学データベースが文献と異なっていたようです。


WATEQ4f で計算していたのですが、文献は PHREEQC でした。扱っている鉱物の数が倍半分ですので、結果も異なったのでしょう。ま、傾向は同じですし、どちらが正解に近いかは、試験室でないとわかりません（局所平衡を仮定し、EQUILIBRIUM_PHASES を使用した Test Run であって、地下水の流れを再現しているわけではありません）。

ZVI は PHASES で Fe を分解し、logK、enthalpy を与えています。ま、その値が分からずに今まで何もしなかったのですが。手元にある参考書には logK しか載っていなかったので、温度変化が扱えなかったんですよね。

とにかく、PHREEQC で再現できたのは収穫です。input ファイルの組み方が確認できましたので、次は PHAST で透過させてみましょうか。

透過反応壁とPHREEQC

以前より、ZVI を用いた PRBｓ の計算をしたかったのですが、情報が少なく、そのままにしていました。

先日より、また気になりだして資料を探していました。

今回、一番に引っかかったのは Amazon。洋書でした。
「なか見!」で本の中を見ましたが、モデリングの詳細は書かれていませんでした。買う必要はないと思い、そこの引用文献だけ検索してみました。（便利な時代になりました）。
引用文献は PDF で入手できました。古い米空軍の PRBs に関する資料です。

これ、当たりでしたね。PHREEQC による簡単な計算結果が載っています。
ただ、肝心な ZVI のモデル化が抜けており、掲載されているinputファイルそのままでは動きません。たぶんこのように組んだはずだと、文献を読みながら作成。で、計算！


結果を見ますと、傾向は合っているのですが、数値が微妙に合いません。

続きは後日。

スタックサイズの変更

並列化済みの Dtransu を動かす場合、スレッド数やスタックサイズは、batファイル内で指定していました。

後輩から、それをソース内で固定できないか？といった問い合わせがありました。
ソース内で固定したことはないのですが、当然可能だろうと触ってみたところ、うまくできません。力不足です。

調べてみると、スレッド数の指定はすぐに引っかかりました。Win7 + Fortran では以下の通り。

call OMP_SET_NUM_THREADS(スレッド数)


しかし、スタックサイズの指定が分かりません。オプションの中を見てもわからなかったので、その日は解決できず。


後日、調べてみましたが、結局明示する方法はわからず。最終的に、Win7側で環境変数を固定することにしました。

OMP_STACKSIZE　を追加すればOKです。スタックサイズはそれほど変更することがないので、固定しておいても問題ないですね。画像では、OMP_NUM_THREADSも追加していますが、ソース内で指定する場合は不要です。

しかし、これらを解決するのに、3夜かかりました。
プロがそばに欲しいですね。

玉石径は3倍？

今日、突然電話がかかってきました。

「玉石の径を推定する場合、ボーリングで確認した長さを3倍するが、これ、何かに書かれていないか？」

前からよく聞かれます。が、使わないので忘れます。せっかくなので、書き残しておきましょう。


玉石径3倍についての記載は、基準類には書かれていません。調べてみましたが、参考書では、以下のp126の一文だけでしょう。
近代図書「ボーリングデータの見方と活用ノウハウ」

HPではいくつかあります。
例えば、東北地質調査業協会でも、3倍の説明があります。
http://www.tohoku-geo.ne.jp/technical/qa/07/index.html

発表要旨もいくつかあります。しかし、3倍で良い、3倍じゃダメ、いろいろです。いずれにしても、発表程度では会検に耐えられません（ココまで指摘されるんですね）。


個人的には、地質によって長短比が異なるため、実際に河原などで測るようにしています。でも、それが3倍よりも小さいと、設計者は3倍として扱うことが多いようです。

ま、そのようなところまで根拠を求める方もいらっしゃると、頭の片隅においておきましょう。

面流量

要素の面を通過する流量について。

境界面なら、節点流量で計算できます。が、モデルの中の節点は+-0なので、集計できません。

計算するとすれば、以下の手順でしょうか？

要素No.チェック
面No.チェック
面の構成節点チェック
節点の座標チェック
3点から２辺のベクトル作成
２辺のベクトルから外積計算
要素流速から合成ベクトル作成
法線ベクトルと合成ベクトルの内積計算
面積をかけて流量とする。


いえ、適当です。組み立てていかないと、正しいかどうか分かりません。書いている途中から実務では無理だと思いましたので。

やはり、ポスト処理のソフトを使うべきですか。

NHKスペシャル 汚染水 ～福島第一原発 危機の真相～

NHKスペシャル　汚染水 ～福島第一原発　危機の真相～
http://www.nhk.or.jp/special/detail/2013/1201/

地下水の流れに汚染がのることは、TEPCOさんなら、重々承知されていたことでしょう。
ただ、わかっていても、できなかったのだと思います。

一方、番組では現在も漏洩の原因を究明していく体制が映されていました。現場で船を操作する映像にも、後ろに余裕のありそうな方が映っていました。

違和感の多い番組でしたね。

浸透流計算が収束しない その2

また、後輩からヘルプ。Dtransu で初期定常が収束しないとのこと。
http://phreeqc.blogspot.jp/2012/08/blog-post_25.html

改良済みの UNSAF なら回るのに、Dtransu だと回らないそうです。地下水や雨がオーバーフローした時の取り扱い方が違うのでしょうね。
Dtransu は、降雨浸透境界に降雨があり、さらに正の水圧がかかると P=0 の水頭固定となります。私は UNSAF を使用しませんので、どのように改良されたか知りません。知る必要がありますね。

原因としては境界条件の張り方や透水係数のオーダーなどいくつかあったのですが、やはり根本は不飽和特性。
これを簡易に回避している資料が、原子力安全基盤機構さんより公開されています。
http://www.jnes.go.jp/content/000117846.pdf

なお不飽和特性に関しては、収束性の観点より飽和度の低い領域では比透水係数の最小値の最小値を0.2として設定して補正を行った(図3.8は補正を行う前のものである)。

なお、AC-UNSAF3DにおいてはDtransu・3D・ELとほぼ同様な結果が得られた。
3D-SEEP及びTOUGH2においては,上述の通り適切な涵養量を設定できずDtransu・3D・EL等にて設定した降雨強度をそのまま使用したため、地表において飽和となった領域に対しても設定した降雨強度を無理やり押し込むと言った計算になってしまい水頭が異常に高い結果となってしまった。特にTOUGH2においては、収束性をよくするための不飽和特性曲線に対する補正が行えず、今回設定した収束条件を満たすまで解を収束させることが出来なかった。

広域解析では比透水係数の最小値を0.2とした補正を施したが、小域は山岳部が存在しないことから負のサクション圧が大きくなることはないと考えられること、メッシュ分割が細かいため解析が収束しやすいこと等を考慮し、小域解析では比透水係数の最小値を0.01に補正した不飽和特性曲線を使用した(図3.9は補正を行う前のものである）。

うーん、0.2は大きいですね。
解決したといえるのかわかりませんが、ま、悩んでいるのは私たちだけではないということです。

有限体積法

有限要素法や差分法はよく聞くのですが、有限体積法は今まで触ったことがありませんでした。

応力など、なめらかな形状が必要な場合にはFEM、地下水などはボクセルの差分で十分だと思っています。FVM はどういった時に使うのでしょうか？周りのソフトは FEM か FDM ですので、昔流行った手法かと思いきや、そうでもないようです。ネットでは、流体解析では FVM というのも見かけます。よくわかりません。

そもそも、有限体積法という名前しか知りませんでした。
調べてみると、以外に容易かもしれません。読んだのは平瀬創也「C#で学ぶ偏微分方程式の数値解法」。
質量保存を考える際、微小体積を仮定しその中の収支を追いますが、FVM もまさにそれでした。格子点周りに微小体積（コントロールボリューム）を想定しすれば、積分形の式を利用しやすくなります。さらにそれを差分形式で近似すればおしまいのようです。
でも、わざわざコントロールボリュームを仮定するのが、まだピンと来ていません。結局差分形式で表現するなら、最初から差分法でも良いような気がします。

非構造格子でも対応できるそうですが、そこまではまだ理解できていません。
そういえば、MODFLOW-USG も驚くような格子になっていましたね。差分法のままなのでしょうか？どうなっているのでしょう。

FVM の解説本は少ないようです。ま、構造格子のみだと上記の本で問題ないので、とりあえずは良しとしましょう。

在宅勤務

今年のはじめだったでしょうか、米Yahoo!の CEO が在宅勤務禁止の意志を示しました。

その後、禁止になったのかどうか知りません。その真意も分かりません。ただ、会社で人と人が顔を突き合わせた方が、ディスカッションが増える、より良いアイデアが生まれるというような意図は理解できます。時代に逆行しているとの論争も出ましたが、個人的には行き過ぎた在宅勤務を引き戻し、モニター越しでなく、面を合わせて話し合うというのは基本的かつ重要なことだとおもます。逆行しているとは思えませんでした。

私の勤めている会社は、在宅勤務を認めていません。

家では仕事をしていないことになっています。

Yahoo!より、進んでいるのかもしれません。

DEM と DTM

国土地理院さんが基盤地図情報を更新されています。

先日のUPで、さらに LP ベースの 5DEM が追加されています。以前、写真測量ベースのDEMしかなかった箇所も整備されていました。いいですね。

そういえば、「DEM と DTM の違いは？」と聞かれたことがります。個人的には、グリッドデータになっているものが DEM という解釈をしていましたが、正確には答えられませんでした。同じ解釈をしている方もいらっしゃいますが、そうでもない解説もあります。
http://gis.stackexchange.com/questions/5701/what-is-the-difference-between-dem-dsm-and-dtm

DTM: bare-earth representation with irregular spaces between points (non-raster).
DEM: gridded raster representation of the DTM.

国土地理院さんの解説では「同義」と解釈してもよさそうです。
http://www1.gsi.go.jp/geowww/Laser_HP/faq07.html

DEMに類似する用語にDTM(Digital Terrain Model)という略語があり、日本語では数値地形モデルといいます。これら2つの言葉は異なる語源からなりますが、同義語として扱われる場合が多く、一般に「DEM」が使われる場合が多いようです。

国土地理院さんの動きを見る限り、取り急ぎ、全国整備が目標なのかと思われます。
今後は時系列で整備していただくと、より便利になりますね。どの地域に何年のLPがあるのか分かれば、崩壊前後のデータがあるのか、すぐに判断できます。
今後に期待です。

河川の流れの基礎方程式

LSFLOW は連続の式（質量保存則）と運動方程式（ナビエ・ストークス）を利用しています。

その詳細は非公開と聞いていたのですが、土研資料にはソースまで載っていますし、砂防学会「地震砂防」という本に、支配方程式の展開や離散化の方法も載っていました。砂防分野では有名なのでしょう。

式の展開を初めて見たときは、独特な定式化だと感じていました。特に、クーロン則などを適用している外力の箇所。解釈は結構、自由なんですね。ま、再現できれば良いわけです。

また、地すべりの運動では鉛直方向の流速が水平方向に比べて無視できるという仮定のもと、式の簡略化が行われています。これを見たときは独特だなあと感じつつ、昔のことですから、計算容量を小さくしたかったのだろうなあ、程度に考えていました。

ところが先日、河村哲也「河川のシミュレーション！」という本を読んでいて、同じような定式化を見ました。河川の流れの基礎方程式というそうです。「鉛直方向は無視」という簡略化は河川砂防分野では常識なんでしょう（私が知らなかっただけのようです）。こういう考え方が河川砂防分野のベースにあったからこそ、LSFLOW も、あのような定式化をしていたのでしょうね。もっといろいろ学ばないといけません。

LSFLOW も学位論文等で FVM に改良されています。LSFLOW3、LSV と呼ばれているコードがそれに該当するのようです。
いずれにしても、inputファイルの作成～計算～可視化（MVS）までの実装は完了しましたので、あとは崩壊前の地形の調整を含めたパラスタのみです。地形の調整は時間がかかるでしょうね。
http://phreeqc.blogspot.jp/2013/10/blog-post_9060.html
もう少し考えないといけません。

Visual Fortran と Win8.1

計算用のPC に Visual Fortran Composer XE 2013 SP1 UP1 を入れています。

先週末、そのPCにて 6コアフルの計算させるていたため、他の作業ができませんでした。良い機会なので、VAIO　＋ Win8.1 に VF を入れて作業しようと思い、インストールに取り掛かることに。
最近知ったのですが、VF のシングル・ユーザー・ライセンスは PC 1台につき１つではなく、1人につき1つなんですね。つまり、個人使用の PC が3台あれば、3台すべてにインストールできるわけです（当然、同時使用はナシ）。

VF をインストールしていると、13/30まで進んだ頃にエラーが発生。何度かやり直してみたのですがエラーで進みません。中途半端に入ってしまったので、uninstall しようとしたら、これもエラー。再起動をかけようとすると、ずーっと再起動中の表示で落ちない。電源ボタン長押しで強制終了し、再度起動する羽目に。
その後、uninstall してから、再度インストールするも、全くダメ。失敗です。
結局、システムの復元でさっぱりしました。

13/30 はちょうどVS 2010 shell をインストールする過程でした。せっかくですので、shell を最新版にすることに。
VS 2013 はexpress 版が3種類ありました。for Windows が2種類。DLページに説明は書かれていませんが、たぶん、ストアアプリ用と、デスクトップ用でしょう。気にせず、for Windows Desktop を DLし、インストール。firewallを落とすのを忘れていたのですが、一発でOKでした。

しかし、VFが2013に対応していないのか、統合できません。
HPで確認すると、Win8、VS 2012までのようです。
http://www.xlsoft.com/jp/products/intel/compilers/fcw/

• オフロードと SIMD 拡張を含む多数の OpenMP* 4.0 の機能をサポート SP1 より追加
• Fortran 2003 規格のサポートを拡張 (ユーザー定義の派生型 I/O) SP1 より追加
• Microsoft* Windows* 8 および Microsoft* Visual Studio* 2012 をサポート

仕方ないのでアンインストールしてから2012を入れることに。2012も express 版を入れてみました。
が、これも統合できません。この時点で朝の4時。続きは翌日に持ち越し。


翌日、気を取り直してPCのメンテから始め、再度 VS 2010 shell 付属の英語版を入れてみました。これはあっさり入りました。こうなると、日本語版も入れてみたくなります。ダメもとで再度チャレンジしてみると、13/30でかなり考えながらも、今度は最後まで到達！何が悪かったのか分かりませんが、入りました。

後々、サポートに聞いてみると、「リリースノートを見なさい」とのお返事。
見てみると、ちゃんと書いてありました。概要は以下の通り。

・Windows XP SP3 ～ 8.1までOK

・IA-3、インテル64 対応アプリケーションのビルドに、Microsoft* Visual Studio* 開発環境あるいはコマンドライン・ツールを使用する場合は、次のいずれか:
o Microsoft* Visual Studio* 2012 2010 Professional Edition 以上、2008 Standard Edition 以降、 2010 Shell、2008 Shell

・IA-32 アーキテクチャー・アプリケーションのビルドに、コマンドライン・ツールの
みを使用する場合は、次のいずれか:
o Microsoft* Visual Studio* Express 2012 for Windows Desktop
o Microsoft* Visual C++* 2010 Express Edition [3]
o Microsoft* Visual C++* 2008 Express Edition 

http://www.xlsoft.com/jp/products/intel/compilers/fcw/2013/Release_Notes_sp1_up1.pdf

2012 PROはもっていませんので、2010 shell で BEST だったようです。そこに行きつくまで、1日かかりました。つくづく、周りにプロがいれば、と感じた週末でした。

Tecplot で水面の座標抽出・差分表示（備忘録）

備忘録です。


① Iso-surface で P = 0m の面（水位）を表示させる。

②“Data” - “Extract” - . “Iso-surface” で水位面を抜き出す（Zoneとして新たに追加される）。すべての zone から一気に抜くことも可。

③“Data” - “Alter” - . “Specify equations”にて、変数を2個追加。名前は何でもOK。ここでは、WL1 （ H を代入）、WL2 （ 0 を代入）を作成。
｛WL1｝={H}
｛WL2｝=0
＊この後、“File” - “Write Data File”で、作成した水位の座標出力可能。
＊ただし、吐き出されるXY座標はZone毎に微妙に異なるので、単純に差分は取れない。

④水位を抽出したZoneでは、それぞれジオメトリが異なっている（Z=H（水位）となっている、XYが微妙に異なる）ため、単純に差分は取れない。そこで、Zone同士を比較する場合、比較先のジオメトリにて、比較元のHを補間してから差分をとる手順になる。まずは補完。
“Data” - “Interpolate” - . “Kriging”にて、Destination zone のWL1に、Source zone のWL1をkrigingで補間・代入。ここでは、Zone146のWL1(=H：水位)をZone74のジオメトリ上のWL1に補間・代入している。

⑤差分計算

Destination zone で指定したzone（ここでは74）を選択し、{WL2}={WL1}-{H} でCompute。

WL1 には Zone146 の値が④の過程で補間・代入されているので、選択した Zone74 の水位との差分が、Zone74 の WL2 に代入される。

⑥Time Stands
差分をとり終わったら、表示したい Zone を Time に追加。

⑦コンター図表示で完成



=========================================================

2013.11.23 追記

kriging 時、Driftはナシです。

Tecplot で水位の差分表示

Tecplot で水面を表示できますので、その差分をソフト内で計算・表示できれば、座標書き出しの必要性が小さくなります。

ま、できるでしょうと思いチャレンジしてみたのですが、できませんでした。水面のみ抽出した2つのZoneを指定し、差分を取るだけなので容易にできそうなのですが、ダメ。いえ、差分を取る計算式を入力する箇所はあるのですが、Zoneを指定すると"incompatible zones "というエラーが出るのです。

ググってみると、フォーラムに以下のような質問がありました。

http://www.tecplottalk.com/viewtopic.php?t=339

I have got a problem with the data -> alter -> specify equations function! I have tried for ages to visualize with this tool an anspect of my flow field, but it doesnt work.
description of my problem:
1. I load two files into tec360 with the same geometry datas but other flow fields. This works, and so i can see 10 zones (5 of each loaded files)
2. What i wanna do now is to subtract the pressure level of the first file of the second one. (zones 1-5 are from the first one, 6-10 from the 2nd one)
3. So i go into the specifiy functions tool and write:
{pdelta} = {p0}[6] - {p0}[1],{p0}[7] - {p0}[2],... (and so on for all the corresponding zones from 1 to 10)
-> but now i cannot compute this value pdelta and visualize it...-> if i try to compute it, there are three error messages which appear:
1. if i don't select any zone in the 'zone to alter' list, the error appears: 'invalid zone set requested'
2. if i select all zones, the following error appears: 'incompatible zones referenced in equation'
3. and when i select only a part of the zones, e.g. for the equation{pdelta} = {p0}[6] - {p0}[1] the zone 1 & 6, the following error appears: 'must alter all zones when creating a new variable'
So my question is, how can i compute and visualize the problem i described? what is wrong with my equation? Is it the syntax or has it to do with my zone selection or status? 

驚きです。
私も、まったく同じことをしようとして、同じところで躓いていました。 欲しい機能もわからないところも同じということです。この質問には解決策が載っていませんでしたので、Tecplot にその機能はないということなのでしょう。マクロを組めば、できるかもしれませんが。

代わりと言っては何ですが、1ステップ毎の計算は可能でした。欲張らないことですね。

http://www.tecplottalk.com/viewtopic.php?t=1124

In order for Tecplot to perform a calculation involving more than one zone (like the time average) all the zones involved have to be "compatible". That is, they have to have the same underlying order, and the same number of data points.

The "incompatible zones referenced in equation" message means that at least one zone had a different oder, or different number or points, or both than the rest of the zones.

For these instructions, think of two incompatible zones, each with five variables (X,Y,Z,T,P), where T is Temperature and P Pressure. And suppose that you want to subract P of zone 1 from P of zone 2:

1) Go to Data/Alter/Specify Equations and create two new variables, by highlight both zones in the Zone to Alter box and entering two equations like these: 

{P1} = 0
{P2} = 0
Then click Compute.
Both of your zones will now have two extra variables. 

2) Still in Data/Alter/Specify Equations, highlight only the first zone in the Zones to Alter box and then erase the two equations you entered moments ago and enter this one:
{P1} = {P}
and click Compute. This will copy P of Zone 1 into P1 of zone 1.
Now you can close out of Data/Alter/Specify Equations.  

3) Go to Data/Interpolate and choose a method. Select Zone 1 as the Source zone and Zone 2 as the Destination zone and P1 as the variable and click Compute. This will interpolate P1 of Zone 1 onto P1 of Zone 2.
P of Zone 1 has now been transferred to 

4) Now close out of Data/Interpolate and go back to Data/Alter/Specify Equations. Erase an equations that are present, and select Zone 2 in the Zones to Alter box.
Then enter an equation like this:
{P2} = {P1} - {P}
and click Compute. P2 of Zone 2 will now hold the result of the calculation of P of zone 1 - P of zone 2

なるほどです。確かに、抽出した水位の座標はそれぞれ異なりますので（水面（Z）だけでなく、なぜかXYも微妙に違う箇所がある）、エラーを吐くわけです。
本来は、上側の質問のように複数のタイムステップ同士を比較したかったのですが、贅沢は言えません。後者で実施してみると、うまくいきました。

続く。

Iso-surface の座標

浸透流計算の結果を Tecplot で表示しています。

地下水位として、圧力水頭 0m の面を表示しているのですが、これの XYZ 座標が抜けません。
施工前後の差分をとりたいのですが、やり方も分かりませんし、座標の書き出し方も分かりません。SoilPlus で計算していたら差分表示も楽なのですが、今回は Dtransu 。 CTC の Dtransu > S+ コンバーターは、私の環境では正常に動作しませんので、可視化ツールに頼るしかありません。
うーん、と悩む時間も勿体無いので、MVS を使用することにしました。

といっても、MVS 用の変換ツールを作るところからです。
データの読み込みは Tecplot 用のソースを再利用できます。書き出し部分だけ、inp ファイルの並びにあわせると、出来上がりです。それでも、バグ取りをしていたら、正常に動作するまで4時間かかりました。

早速、書き出した inp ファイルを MVS へ読み込んでみました。
が、表示できません。汚染の表示のように、Iso-surfaces として圧力水頭分布を表示できると思ったのですが、上手く行きません。こちらも、うーん。変換をまだ誤っているのかもしれません。ま、Iso-surface を表示できても、座標書き出しまでたどり着くには、まだ考えないといけませんし。こちらも、先は長そうです。


食事休憩をとった後に Tecplot に戻ってみると、なぜかあっさり Iso-surface の座標を書き出せてしまいました。なんだったのでしょうか。後日、備忘録をつけておきましょう。

ま、こんな日もあるでしょう。

三軸圧縮試験の拘束圧 その2

試験機に詳しい先輩に聞いたところ、使用しているロードセルの精度は1/4000。最大荷重は2kN。この値だけであれば、低拘束圧でも大きな問題は生じないように思えます。
が、先輩が言うには、低い荷重領域はやや精度が劣るとのこと。

例えば、直径5cmの供試体に、400N（最大荷重の20%）をかけると、約200kN/m2。低拘束圧下での破壊を考えると、やや大きな値です。背圧を100kN/m2かけても、クリアできない試料もあると思います。
では、背圧を200kN/m2にすれば良いのでは？と聞くと、それも経験上ダメとのこと。理論上はOKでも、経験上、静水圧程度にしておかないと、（実際は違うかもしれませんが）構造が壊れやすい？そうです。

他社のプロとも話しました。試験室でロードセルとPC画面を見ながら話をしたのですが、こちらも似たような精度でした。どちらかというと、不均質性の影響がより大きくなるだろうと。ま、その通りです。サンプリング・運搬・トリミングによる乱れ・供試体の不均一性はそのまま出てくるでしょうね。室内にだけ精度にこだわっても、ダメということです。

御二方と話して得た結論としては、「低拘束圧のみで試験を実施するなら、ロードセルを変えたほうが良い」ということ。うーん。ま、そこまでこだわれということかもしれません。


あと、先輩、私、プロ三者に共通した意見がありました。
「低拘束圧のほうが、cが出るだろう」
基準では逆の説明がされています。粒子破砕でも起これば、高拘束圧のほうが寝てくると思いますが。どうなんでしょうね。

試験一つとっても、まだまだ分からないことが多い状態。一歩づつ、片づけていきましょう。

三軸圧縮試験の拘束圧

三軸試験の拘束圧について後輩から質問がありました。

三軸試験の拘束圧はc・φを決めるのに重要な要素です。が、多くの方は気を使われていないようです。現場担当・試験担当などの分業化、試験の基準化の弊害でしょう。現場が「三軸試験をお願い」と言えば、基準に従った結果は上がってきます。が、側圧の詳細な設定方法は基準に載っていませんので、試験者によってはゲージの読みやすい値や習慣で、とんでもない拘束圧を設定することもあるでしょう。拘束圧の設定方法に疑問を持った（設計者）は「いいね!」です。

余談ですが、他社の熟年設計者が、お客様に対し、拘束圧をどの程度に設定するか聞いているのを見たことがあります。お金をもらってコンサルティングする立場のはずですが。これは逆に論外です。

拘束圧に関して、基準に載っているのは下記ぐらいでしょうか。

国土技術研究センター「河川堤防構造検討の手引き」H24改訂版

http://www.jice.or.jp/siryo/t1/201202090.html

従来、三軸圧縮試験および一面せん断試験における拘束応力の設定に配慮不足な面があった。すべり面計算に用いる三軸圧縮試験等のせん断強度試験は、発生すると予想されるすべり面の深さにおいて、発揮する強度が評価されるように、低い拘束応力範囲を含むように設定することが望ましい。高拘束圧下の試験結果から得られた粘着力を見込むと、低拘束圧下で過大な強度となり、過大な安全率が得られることもあるため、粘着力の評価に必要な注意事項である。

ただ、具体的な設定方法については書かれていません。明確に書かれた基準は現段階でないと思います。

本来、盛るのか、掘削するのかによっても変化しますし、安定計算時のすべり線が決まっているかどうかによっても変わります。扱うモデルが2次元か3次元か、欲しい値が過圧密の範囲内なのか正規圧密状態なのかでも異なります。現状ではそれらを書き出すのが面倒なため、常識？の範疇となっているのでしょう。

先日、他部署の先輩と話していた時に、ある資料の存在を知りました。またも河川、下記土研の解説書です。この分野、河川は強いですね。
比較的初心者向けの解説書と思われますが、その中に拘束圧について書かれている箇所がありました。公的な書き物で具体的な記載があるのは初めてではないでしょうか？
土木研究所　地質・地盤研　土質・振動チーム「河川堤防の浸透に対する照査・設計のポイント」
http://www.pwri.go.jp/team/smd/topics-seepagepoints.htm

ただし、これにも問題があります。
最大で50kN/m2だと、他は10, 30kN/m2程度でしょうか。拘束圧は設定できると思いますが、10kN/m2はかけたことがありません。最低でも20kN/m2程度です。機械にもよると思いますが、軸方向の精度が保てるかどうか心配です。
あと、モール円が詰んでしまい破壊線を引き難く（c・φを決定し難く）なりますね。ま、これは(有効)応力経路のグラフにて、破壊点に対し最小二乗近似を取ればクリアーできますが。

続きは後日。

崩壊の原因は？

今年は台風が多いせいか、土砂災害も多く発生しました。

今年、呼び出されて見た中で、印象に残っているのが、ある切土の崩壊。
岩盤の上に厚い扇状地性堆積物が載っていたのですが、そこをオープンカットで通過していました。で、台風により上位の堆積物だけが崩壊。原因は岩盤との境界を流れる地下水です。土工指針や教科書に載っている、典型的な崩壊でした。

先日、お世話になった方から執筆図書の案内があり、ネットで注文しました。届いた本を見ていると、学生向けの内容でした。が、そこにも典型的な崩壊例として、岩盤上の土砂の崩壊が載っていました。なぜ、設計者がそんな初歩的なミスを犯したのでしょうか。

ま、少しわかる気がします。たぶん、仕事として淡々とこなしてしまったのでしょうね。

以前、掲示板（雑談板）で、高級レストランか、町の定食屋化か、どちらのような会社を目指すべきかという話題を掲げた方がいらっしゃいました。詳細は忘れましたが、高度かつ独自の技術で勝負する会社か、小さな問題でもすぐに対応できる声のかけやすい会社か、ということだったと思います。
私は、どちらも目指すべきだと思います。逆に、心が入っていれば、どちらでも良いとも思います。どんな仕事も人と人のやりとりです。心がなければいけません。今話題の食材誤表示（詐欺）と同じくらい、インフラにも問題は潜んでいるかもしれません。

確率降水量

後輩二人から、確率降水量の設定についてレクチャーを受けました。

雨量の引伸ばしを自分で行うのは初めて。ですが、確率降水量の考え方や算出方法自体は、他の方がされているのを今まで見てきたので、ま、できるでしょう程度に考えていました。降雨データを用意して、確率分布関数（確率密度関数）を複数あてはめ、最も適合したものを選ぶというだけです。気象庁のHPが分かりやすいと思います。
http://www.data.kishou.go.jp/climate/riskmap/cal_qt.html

ですが、できませんでした。
どのような雨量データを用意すべきか、という初期の段階で躓いたのです。
河川の方なら台風など一雨で期間を区切ると思います。が、私の適用したい現場では、数か月～1年のオーダーで雨量の引伸ばしを行う必要がありました。そこで後輩の出番。1時間ほど話ながら考えてもらい、最終的には解決しました。ありがたいことです。

ちなみに、計算に使ったソフトはコチラ↓
一般財団法人国土技術研究センター「水文統計ユーティティ」
http://www.jice.or.jp/sim/t1/200608150.html

分布関数は16種類。計算はほぼ自動。若干、使い方に戸惑うところはありましたが、よくできていると思います。こういったソフトを無償公開されるのはありがたいことです。今後も制約なく使い続けることができれば良いですね。

ビジネスモデル特許

ビジネスモデル特許というのは、結構身近にあるようです。

携わっている仕事の一部が、どうやら他社のビジネスモデル特許に引っかかりそうだと分かりました。いえ、先輩はずいぶん前に御存知で、その辺りのことは既に手を回し、処理図済みだったのですから驚きも倍です。技術以外のことも勉強しないと、やっていけないようです。

今回は、例えば以下のようなことでした。

ある先生が、FEMをつかった変位解析手法を論文で発表しました。それを見たA社がコーディングし、誰でも使えるように公開しました。B社はそれを改良し、変位についてより精度を高めることに成功しました。これも論文で公開され、皆が自由に手にすることができ、技術の発展に大きく寄与しました。
ここで登場したのがC社。B社の改良ソフトをある現場条件に利用し、変位ではなく応力についての利用法をビジネスモデルとして特許申請しました。今まで、皆が変位にばかり注目し、たまたま一緒に求まる応力について触れていなかったのです（実際はそんなことはないでしょうが）。
それを知らないD社が、B社の改良ソフトを使って応力について求める業務を行いました。が、運悪く、その条件がC社の申請した条件と合致してしまいました。

似たことことが、現実にあるんですね。
誰が悪いのかというと、感情論ではC社です。が、違うんでしょうねえ。

そういえば、IPS細胞の作成等でノーベル賞を受賞した山中教授は、特許管理についても力を入れていると、テレビで見たことがあります。知財管理・契約を専門とする部署があるようで、特許の権利を大学に取得・管理させることにより、学術研究への無償提供による発展、他者や企業の特許取得による実用化の遅れ・独占防止、研究分野全体への支障を防ぐというような趣旨だったと思います。

ま、これは本当の特許で、ビジネスモデル特許とはかけ離れたものなのですが、通じる物はあると思います。研究者や企業は良い論文を発表するだけでなく、その発展や実践も視野に入れ、先々手を打っておく必要があるのでしょう。

不合格通知

H25 技術士試験の不合格通知が届きました。

選択ＩＩがBでした。理由は明らか。選択問題が少なくなり、実質、選択の余地がなくなったこと、そのうち1問が60%以上回答できるレベルに達しておらず、専門外の問題に手を出したこと。
http://phreeqc.blogspot.jp/2013/08/blog-post.html

他人から「お前はダメだ」と言われると、客観的にダメだとよく認識できます。この年になると、なかなか面と向かって「ダメ」と言われませんので、これだけでも良い収穫でしょう。

先輩に、専門外の問題を答えたことを話すと、「口答試験で、その専門外の方が来たら最悪」と言われました。なるほど。
受験以降、試験対策として専門外も勉強する必要があると考えていたのですが、それでは口答試験の時にタジタジになる可能性が出てきます。やはり、専門分野に絞って勉強する方が良いのでしょう。

会社の掲示板に合格者名簿が出ていましたが、今年は10名強と非常に少なかったですね。異例です。試験制度の変更のため受験を控えた方が多いことに加え、社内の高齢化、金銭面での社内補助の減が要因と思います。残念ですね（という私も久々の参加でしたが）。

ま、60点でOKの試験です。それに満たないようでは、その分野の技術者と名乗れません。まだ見習いです。せめて技術者としてのスタートラインには立ちたいものです。
来年も短期目標としましょう。

透水係数と回転マトリックス

3次元の浸透流解析でパラスタをしていました。

なかなか合わない個所があり、kx, ky, kz を変更してみることにしました。それが正しいかどうかは別として。

砂岩・泥岩の互層であり、山を歩いた感じでは一様に傾斜しています。その平均的な走向傾斜は把握しています。
また、Dtransuであれば方向余弦の入力で回転した透水係数(軸)も計算に取り込めます。マニュアルに記載されていますし、2年前に確認しています（なぜ符号をそのようにしたのかは全く覚えていません。GTRANの座標系だったかな？）。
http://phreeqc.blogspot.jp/2011/07/dtransu.html

で、いざ反映しようと思うと、あれ、分かりません。方向余弦を出すにも、2軸で回転させてしまうと、わからない角度が出てくるのです。回転マトリックスを使うにしても、平面上を1回転までは計算できるのですが、次に傾斜方向へもう一度回転となると、その方法が分かりません。どちらかの軸だけで1回回転ならOKですが。

調べてみると、走向傾斜はオイラー角で表現するのが便利だとわかりました。これ、知りませんでした。回転後の軸を使って、もう一度回転させる、そう考えれば回転マトリックスを2回使うだけ、積をとるだけでOKでした（数学、もっとやっておくべきですね）。個人的には、こちらのサイトが分かりやすかったです。
http://www6.ocn.ne.jp/~simuphys/daen1-1.html
wikiのGIFがイメージしやすいです。
http://ja.wikipedia.org/wiki/%E3%82%AA%E3%82%A4%E3%83%A9%E3%83%BC%E8%A7%92

上のサイトでは3回転させてますが、走向傾斜であれば2回転でOKなので、もっと簡単な式になります。

この程度ならEXCELで配列数式があるかと思いましたが、ないですね。ま、それほどややこしいものでもないので、sin、cos と手入力で仕上げました。走向傾斜を入れたら、方向余弦を算出するようにしました。もっと簡単な方法もあると思いますが、ま、今回はこれでOKです。


で、目的だった走向傾斜と透水係数の異方性を考慮した計算をしてみました。

結果、今度はほかの箇所が合わなくなりました。うーん。



*********************************************************
2013.11.02
以下にEXCELデータをUPしました。
https://sites.google.com/site/geochemist001/resources/RMatrix

計算コードの誤り

先日の、どこが改良されたのかわからない Fortran のコードを読んでいました。
http://phreeqc.blogspot.jp/2013/10/avsmicroavsmvs.html

説明書もコメントもありませんので、数人の書いたコードを順に比較するしかありません。

最初の内は、軽微な修正で確認も早く終わりました。が、最後の方で、いくつか重要な修正がありました。変数の追加や配置、計算内容を見ると、何をしようとしているのか意図はわかります。が、どうも式がおかしい。

半日くらい考えて、たぶん間違えられたのだろうと、と結論付けかかったころ、コードが入っているフォルダの中にReadMe.txtを見つけました。お、っと思って開いてみると、ビンゴ。変更点の概略と参考文献が書かれてました。


文献を確認したが、やはりコードが違っていました。危ういですね。

といっても、その方が書かれた前の段階で、既に誤っている箇所がありました。計算部分ではなく吐き出しの箇所なので、修正は簡単です。
私の使用している post で読みやすいように、あわせて修正しておきましょう。

地震関連の調査

先日、地元の方に事業説明をしていたのですが、事前の情報が誤って伝わっており、「活断層の調査かと思っていた」と言われました。近いものはあるのですが、活断層調査ではありません。

そういえば、以前お世話になっていた方に原発の調査 = 活断層調査に声をかけられました。が、無理。おもしろそうだったんですけど。


全国で地震・耐震がブームになっているようです。

ブームで終わらせないたくないですよね。
点検や補強、数十年耐える物を作るのも良いですが、それだと一過性になりそうです。ヒトは必ず忘れます。できれば子供達への教育にお金を回したいですよね。
ベースアップした技術、というよりも意識で、さらに次の世代へつないでくれることになるでしょう。

崩壊前の地形

時間があったので、崩壊後、土砂がどの程度移動するかを再現していました（今度は LSFLOW です）。

最初は簡単に考えていました。

崩壊前後の LP を使用すれば、モデルは簡単にできてしまう、と考えていたのです。が、これが誤り。多くの場合、崩壊後の LP はありますが、崩壊前はありません。写真から復元したとしても、その精度はLPに比べて落ちます。場合によっては10m以上の標高差が生じます。今回も部分的に5～10m、最大20m近いギャップが生じました。

http://phreeqc.blogspot.jp/2012/11/5m_4.html

http://phreeqc.blogspot.jp/2012/11/5m-2.html

崩壊後のLPにしても、そのままではすべり面として使えません。冠頭部はそのままで使えるのですが、末端部は土砂で隠れています。そのため、ボーリング結果を使って面を作ってやらないといけません。こんな当たり前のことも、手を動かすまで気づきませんでした。迂闊でしたね。

こういった前処理に手間をかけることが、精度を大きく向上させるための一因になります。

案外、パラスタより時間がかかるかもしれません。

動的配列

EXCEL2010 の VBA で、TXT の並び替えをしようとしました。

あるソフトの吐き出した結果が、以下のような並び。
X1 Y1 Z11
X2 Y1 Z21
・　・　・

それを、並べなおして、新しい TXT ファイルにしたかったのです。
1 1 Z11
1 2 Z12
・　・　・

TXT ファイルを開いて、ソートして（他にもいくつか計算過程を入れて）、書き出すという流れにしようかと思ったのですが、X1, Y1 を1から始まる番号で書き直すので、Zのみの2次元配列でこなした方が早いことに気づきました。で、以下のように記述。

Dim Z(1 To NX, 1 To NY) As Long

ところが、「定数式が必要です」とエラーが出ます。
調べてみると、通常の2次元配列に2変数が使えない仕様のようです。知りませんでした。
動的配列の場合は、以下の通り書けばよいようです。

ReDim Z(1 To NX, 1 To NY) As Long


微妙です。が、忘れそうなので書き残しておきましょう。

物理演算の手抜き

3ds MAX を触っていた時に、NVIDIA の PhysX が使えるということを知りました。

MAX で PhysX の何が使えるのか、精度はどうなのか、など、知りたいことは多いのですが、なかなか詳しい解説に出会えません。が、GPU を計算に使用でき、レンダリングも可能な環境が既に用意されているのなら、使わない手はありません。少しづつ調べましょう。

先日の PBF は PhysX に取り込まれるかもしれません。
http://physxinfo.com/news/11109/introduction-to-position-based-fluids/
http://phreeqc.blogspot.jp/2013/09/blog-post.html

粒子法の SPH は既に実装済みのようです。でも、MAX でこの SPH が使えるのかは知りません。
http://sa08.idav.ucdavis.edu/CUDA_physx_fluids.Harris.pdf

ただ、引っかかっているのは「手抜き」が入っているのかどうか。以前見た本で書かれていました。
http://phreeqc.blogspot.jp/2010/11/cg.html

先日見た工学社「物理エンジンPhysX&DirectX10」にも、以下のように書かれていました。

「物理エンジン」は、与えられた時間内に処理を完了するため、必要な計算を捨てることがあります。

物理計算を”手抜き”します。

「物理エンジン」の物理学は、不完全な場合が多くあります。

でも、どのエンジンが、計算中のどこを、どの程度「手抜き」するのかが分かりません。
こういうのは、CG分野で学ばないとわからないのでしょうか？

個人的には、物理法則に忠実に問題を解いた方が better だとは思います。が、扱っている多くの問題では、所詮、パラスタで再現しているだけ = パラメーターの物理的根拠がそれらしければ良い = 同定値が真理かどうかは別問題、なので、多少手抜きしていようが再現できるのであれば良いのではないかと思います。
ま、話す場所と、相手は選ばないといけませんが。

こういうのを読むと、手を出そうかどうか、迷いますね。

微分形と積分形 その2

図書館で流体力学の本をいくつか借りてきました。

今回、役に立ったのは以下の本です。

松本洋一郎監修・山口浩樹著「道具としての流体力学」

基本的な式の展開が丁寧に掲載されています。どちらかというと、専門書というよりは一般向けかもしれません。そのため、全体的にやや物足りない感はありましたが、知識の欠落していた部分は簡単に補うことがでいました。

特に、積分形の併記が良かったですね。
質量保存則や運動量保存則より PDE を導く際、私は今まで微小領域の収支より展開する形しか知りませんでした。この本では、積分形も併記されており、2つの方法があることを知りました。これはラッキーでしたね。最近読んでいた論文が、この積分形の展開を使用していたのですが、なぜこのような数学的な展開をするのだろうと思っていたのです。これも一般的なんですね。基礎が抜けているというのは、恐ろしいことです。
微分形は微小領域（直方体）の収支から、テーラー展開経由で PDE を構成して行くのに対し、積分形は直方体にこだわらない領域の質量や運動量を体積分や面積分を使って表現し、展開して行きます（最終的には微分形で同じ式にできます）。微分形・積分形というのは、単純に式の形だけでなく、このような2種の展開の前提をも指すのかもしれません。


ナビエ-ストークスの解釈も問題ないようでした。導出も簡単です。

これで、ようやくスタートラインに立てました。

サラリーマン技術者の報酬

私の長期出張が終わったのと、次の仕事を一緒にする方が支社に来られていたのとを兼ねて、数人で飲んでいました。

いろいろ話をしていたのですが、驚いたのは報酬の話。

つい先日、私の基本給が2年間で6千円しか上がっていないことに気づきました。ここ2年間、ノルマの1.3～2倍の成果を上げ、利益を数千万出していたにも関わらずです。驚きでした。10年前と比べても、下がっていました。

これを肴に飲んでいますと、「それはまだ良い方」とのこと。　？とおもって聞いていると、皆、千円上がっただけだったり、下がったりとのことでした。世知辛い話です。

今年度の12月末までの作業量を数えると、ざっと100人。一緒に飲んでいる方との次の仕事を入れても、心身ともに余裕でこなせる量です。しかし、許可された残業＋就業時間ベースで考えると、80人程度しか費やせません。数字上は業務過多となっています。さらに先日、管理職より追加の仕事を打診されました。が、断りました。受けたとしても、サラリーマン技術者としての報酬は0円です。管理職による部下の管理もこの程度ということです。（ま、ここ2年の評価（報酬）に気づいていなければ、受けていたでしょうが。）

先日、NHK でブラック企業の話題が出ていました。

http://www.nhk.or.jp/gendai/yotei/index_yotei_3403.html

人の世です。それほどきれいな世ではありません。

F の構成

圧力項は、浸透流の動水勾配の項と似ています。

勾配が正なら逆方向に水が流れます。離れた位置の水頭が高く、動水勾配を正とすれば、逆方向に水が流れてくるので負。力も負。だから勾配にマイナスがついていると解釈すれば良いでしょう。速度勾配が大きい（正）ほど、逆向き（負）の流れが強くなるということを示している項です。

粘性項は先日の通り。プラスということは、抵抗ではなく、速度差によって微小体積に発生する力ということでしょう。同じ速度なら、ネバネバしている方が力を受けやすいというのは感覚的にOKです。
違いますね。粘性項は速度勾配の微分（速度の2階微分）なので、プラス方向の速度の差（右に流れるのを正とすれば、流速分布が右に凸な状態）が大きければ大きいほど、負になります。つまり、言葉の通り、内部の粘性抵抗ですね。

外力はプラスで問題ないでしょう。重力がZ軸の負に働けば、負の方向に速度が大きくなります。


これでナビエ-ストークスの右辺に示された3項について、特に疑問は発生しなくなったわけですが（といっても正解を見たわけではありませんが）、この3項だけで F が構成されるというのがよくわかりません。すっきりしませんが、ま、流体の多くはこれで説明できるので、そういうものだうことでしょう。

これで、SPHの前提も、今読んでいる論文の前提もイメージできました。

ラグランジュ微分のオイラー記述

∂u/∂t + (u・∇）u

ナビエ-ストークスの左辺は上記のように書かれています。これについては下記のサイトがよく整理されていると思います。
http://sky.geocities.jp/ima_ich/

Lagrange 微分の Euler 的表現
http://www2.kobe-u.ac.jp/~iwayama/teach/gfd/2012/chap2.pdf


速度のラグランジュ微分 Du/Dt （加速度) をオイラー（空間）記述しているので、移流項が出ているということでしょうねえ。この辺は、数年前に習いましたので、たぶん、あっているでしょう。

昨日の粘性項と合わせれば、移流拡散と同じ型です。スカラーかベクトルかの違いはありますが（大きな違いですが）。さらに、ナビエさんは圧力項がついているというだけのようですね。

そうすると、ナビエ-ストークスの式全体としては、加速度=力の関係になっています。右辺に1/ρがついていますが（粘性項は、動同粘性係数=粘性係数/ρとして消えていますが）、これを左辺にかけ戻せば、密度×単位体積にかかる加速度と考えられますので、ρ×V = m で ma=F、つまりニュートンの運動方程式を記述しているということでしょう。

この F として、粘性、圧力、外力を考えているということですか。どうやらここがポイントのようですね。

粘性項

ナビエ・ストークスの式で、まず引っかかったのが粘性項。

ν∆u

数式としては移流拡散の拡散項と同じ形ですし、スカラーがベクトルになっているだけです。ちなみに、2階編微分の出所は、拡散と同じで微小領域の差を取っているだけでしょう。個人的には、Dtransuのマニュアルの絵が思い出しやすいですね。

導出としては、たぶんこうだろうと、予想がつくのですが、イメージができません。拡散量はxi方向の濃度勾配に比例定数Dがかけられたもの、粘性の場合はxi方向の速度ujの勾配に比例定数がかけられたもの。うーん。どういうことでしょう。せん断方向の速度差が大きくなるほど大きくなる力が粘性ということでしょうか？

ググってみると、わかりやすい説明がありました。感謝です。
http://hb3.seikyou.ne.jp/home/E-Yama/Fluid7.pdf

運動している流体で互いに接している層の間にズレが起これば接線応力（せん断力）が現れるような流体を粘性流体と呼びました。 

ニュートンの仮定壁面よりy, y + dy の距離における速度をそれぞれu, u + du とするとy + dy にある層がy にある層をx軸方向に引っぱる力F は，この2 層の接触面積A とdu に比例し，dy に逆比例する。式で表すと比例係数をμ として

F / A=μdu/dy

これらの仮定はいろいろな実験結果より正しいことが認められています。この仮定に従う流体をニュートン流体と呼びました。
右辺F/A は，単位面積当たりの力でこれをτ で表すと

τ = =μdu/dy

となり，τ をせん断応力（Sheering Stress），比例係数μ を粘性係数（Coefficient of viscosity）と呼んでいます。また，粘性係数μ を密度ρ で割ったものを

ν=μ/ρ

で表し，ν を動粘性係数（Kinematic viscosity）と呼んでいます。

なるほど。定義なわけですね。
引っ張る速さの差が大きくなるほど大きくなるのがせん断力（抵抗力）、せん断力が現れるような流体が粘性流体、それが線形で成り立つ流体がニュートン流体、その係数が粘性係数ということでしょう。いたってシンプル。
ハチミツのように、ネバネバした液体をスプーンで早くすくうときに受ける抵抗力のようなものでしょうか？早くすくうほど、受ける抵抗が大きくなりますので。水の中で動く場合に受ける抵抗も同じですね。飛行機が飛ぶときに受ける抵抗も？要は、ネバネバして粘性係数の大きな流体や、流体中でなんらかの速度差が大きく生じる場合に取り扱うべき項と解釈しました。（流体力学の本を読まないと正解はわかりませんが）

というか、これ、昔、習っていますね。完全に忘れていますが。

ま、イメージはできたので、良しとしましょう。

微分形と積分形

式を追っていくと、「これを微分形で表すと」という表記が出てきました。

がっつり体積分の式だったのですが、いきなりそれを微分形に直すというのは意味が分かりません。過程も出ておらず、結果しか示されていませんので、変換のヒントすらありません。
そもそも、私が不勉強のため、「微分形」という用語（とその変換）があることすら知りませんでした。もっと、若いころに物理や数学を学んでおけばよかったと、改めて後悔。

で、ググってみますと（出張中なので参考書はありません）、比較的よく使われる用語のようでした。見たところ、右辺・左辺ともに体積分で統一したら、中身は一緒でしょ、その中だけ取り出せば、∇があるので微分形です、というようなことだと勝手に理解しました。ま、∇でも∂でも微分があれば微分形と呼ぶのでしょうねえ。真の定義はわかりませんが、直面している式は、その手で解決できそうです。

左辺が偏微分のついた体積分、右辺がF・nを含む面積分です。右辺をガウスの発散定理を使って体積分に上げてやると、∇・Fが生まれて、中身が一緒でしょ、はい、左辺に∂/∂tのつく微分形のできあがり。ということで、展開の結果は合いました（用語の使い方はよくわからないので、帰ったら調べないといけませんが）。

次数を上げて解きにくくしているように見えるのに、微分形の簡単な式になる、という展開の発想がすごいですね。面白い。また、物理的な意味は同じで微分形・積分形の両方で表現できるというのも（私にとっては）新鮮。

ナビエ・ストークスまでは、もう少し時間がかかりそうですが、ゆっくりいきましょう。

ガウスの発散定理のイメージ

論文を読んでいると、ナビエ・ストークスの式が出てきました。

恥ずかしながら、今までナビエ・ストークスを使ったことがなかったので、導出を見てみようと思い調べ始めました。

が、最初に出てきた質量保存式の意味が分かりません。いえ、今まで成分を書き下した式を見慣れていたので、ドット・プロダクトを使った表示の物理的意味がスッと入ってきませんでした。で、逆に成分を書き下して、体積をかけて、流速（流量）の差を取る式に直して・・・・と、イメージできるまでさかのぼって確認しました。やはり使わないと忘れますね。

以下も忘れていました。

発散 (div)   ：∇をベクトルに、形式的に内積として作用させる。∇・u
勾配 (grad) ：∇をスカラーに、形式的にスカラー積として作用させる。∇φ
ガウスの発散定理：∫∫∫v ∇・u dV = ∫∫s u・n dS

質量保存式で流入量と流出量の差を整理すると、∇・u （= ∂ui/∂xi） に行き着く。これは微小体積の発散（湧き出し量）。ガウスの発散定理の左辺は、それを（求めたい領域で）体積積分したもの。u を法線方向に直すと、微小面積当たりの流出量。それを領域の表面で面積分して湧き出し量にしたのが右辺。当然、左辺と右辺は同じ。大雑把にはこんなイメージで覚えておくと忘れにくいでしょうか？

忘れるでしょうね。若いころに叩き込んでおくべきでした。
ナビエ・ストークスは後日です。

AVS MicroAVS MVS

古い土砂移動のコードを動かしていました（LSFLOWではありません）。

数人が改良を重ねているようなのですが、どれが最新版で、どこを改良したのかコメントがありません。ただでさえコメントのないソースは、イマイチよくわからないのですが、とりあえず dimension や変数を拡張し、大きなモデルでも動くようにしました。


結果がAVSで読める形で吐き出されます。
しかし、タイムステップ毎に fld データとして吐き出されているため、アニメーションにするには一手間必要です。

まず思ったのが、AVSの利用。
ですが、他部署の方が休みのため、使用できず。

次に考えたのが、MicroAVS。
最新版がないので、試用版を落とし試しました。
補助ツールとエディターで、時系列のfldを作ってしまえばOKです。これは簡単に動画まで持っていけました。

次、MVS。
AVS の拡張が EVS, MVS なので、吐き出された fld が load EVS field　モジュールでそのまま読めるだろうと。しかも、モジュールを組んでおけば、MicroAVSより手間がかからないのではないか、と思い試すことに。
しかし、吐き出された個々の fld データや、タイムステップを記述した fld データは読めるのですが、アニメーションにできません。TCFを作っても、ステップを指定すると落ちてしまいます。当然、MicroAVS の v ファイルは読めません。ヘルプを見ても、fld を順番に読み込んで動かす方法がわかりませんでした。

うーん。ま、今日はここまで。
とりあえずはソースを動かし、アニメーションにできたので、ゴールにはたどり着けそうです。

山の粘土

離島での調査をしていました。

今回は山の調査なのですが、なぜか海成粘土のような盛土が出現。海の粘土を締固めて盛るというのは聞いたこともなく、しかし、山の粘土でもなさそう。

迷っていたので研究職に相談したところ、「分析よりも、残渣を見ろ」とのこと。

確かに。コア観察をサボってはいけません。


で、蒸発皿を現場に持っていき、コアを取って水洗い。ルーペで観察を続けましたが、貝殻片や黄鉄鉱は入っていません。結論から言うと、海のものではなく、淡水のもの。

しかし、どこから、どのように堆積した粘土を持ってきたのでしょうか？結局、そこまではわかりませんでした。


常日頃、1m下の地盤や岩盤に振り回されているレベルです。まだまだです。

物理演算と粒子法

3ds MAX の massfx, mparticles を触ってみました。

触ってみて初めて気づいたのですが、CGの分野ではアニメーションを作るのに、手で（指示で）ものを動かすのですね。そりゃ、物理演算が可能であれば、そちらの方が容易にアニメーションを作成できるのでしょう。やはり、この分野のほうが、進むはずです。

以下の理論も、途中から理解できなくなりました。動画は今年4月ごろに話題になっていたようです。別の動画には GTX 580 レベル（張り付けた動画は680のようです）でリアルタイム演算ともありましたが、すごいですね。
http://mmacklin.com/pbf_sig_preprint.pdf

下側は Rigid Body との組み合わせですが、土の物性を与えた粒子に置き換えておけば、土砂の流れもシミュできそうです。

土石流や深層崩壊の場合、最初は土や岩石で、崩壊後はほぐれた土砂が流体のように移動し、停止すれば土（堆積・天然ダム）になり、水は流れ去る、といった mass としての物性の変化（各相の物性自体は変化していません）があります。私は土石流や天然ダムのプロではないので知らないだけなのかもしれませんが、それらを1相として計算する場合、土砂の崩壊や落石の集合体に近いイメージで取り扱っているのでしょうか？

mass として簡略化するよりは、固相に水を与えて崩壊させるなど忠実にモデル化する方が、結果的には近道なのかもしれません。

ま、私には遠いゴールです。少しづつ、追っていきましょう。

粒子法 その2

最近、また粒子法について、手を出そうかと思い始めました。

2年前、概要は把握しましたが、それ以降は手つかずの状態でした。
http://phreeqc.blogspot.jp/2011/08/blog-post_12.html

で、再度、教科書を読んでいるのですが、やはりそれ以上はダメ。
SPH の基礎理論はわかります。数式も追えます。MPS は概要のみ。数式の一部がまだ理解できていません。離散化はOK。でも、それを実装するのは全くダメ（CPUのみなら、FEMより簡単だと思うのですが）。そんな状態です。要は、理解していないということです。

これ、ゴールの遠いところが、腰の重い原因の一つでしょうね。
計算以外にも、GPGPU やメニーコア対応は当然の分野。ポスト処理としての陰関数曲面の作成や、レイトレーシングなど、数値計算とは異なった知識も必要です。手に入る情報は簡単か難しいの両極端。市販ソフトも売られている現状、それらを一つづつ追う気力が出てきません。3ds MAX などの物理演算で何とかならないか、と、近道を探したくもなります。

ま、いずれも言い訳なのですが。


土木分野で粒子法がメジャーになるのは、幸か不幸かまだ先でしょう。良い手法ですが、変化の遅い分野ですので。
焦らずに、ゆっくり実装していきましょう。

3軸圧縮試験 供試体の直径

一軸や三軸の供試体の直径は35mmか50mmが多いと思います。

地盤工学会基準はφ35～100mmですが、φ35mmやφ50mmなら乱れの少ない試料を成形し、標準的な三軸圧縮試験機で対応できます。

先日、先の技術者と話をしていて知ったのですが、中型三軸圧縮試験機では、φ100mm以外にφ75mm（会社によってはφ65mmも）のペデスタルがあるそうです。基準に、以下の記載があるので、シンウォールやトリプルの内径に合わせたペデスタルが作られているとのことでした。

「側面の成形によって過度な乱れを与えるおそれがある場合は，サンプルチューブから取り出した試料の側面の成形を省略してもよい。」

教えてくださった方は「最近の流行」と言われていましたが、気を使ってくださったのでしょう。φ75mmは結構有名なようですので、昔からあったのでしょうね。そりゃそうです。トリプルで採取した試料を、わざわざリモールドでφ100mmにすることはないですから。

乱れの少ない試料として採取できるのが基本は砂質土まで。特殊なサンプラーを使わない限り、トリプルで礫質土の採取は困難です。したがって、最大粒径10mm程度までは対応可能なφ50mmの標準的な試験機で、通常は事が足りるといった状況になります。そのためか、会社には中型三軸圧縮試験機がないので、φ75mmのペデスタルは使ったことも見たこともありませんでした。

先日、試験を見せていただいたのですが、やはり大きいですね。供試体の高さが直径の2倍以上ですから、φ75mmで最低15cmになります。4供試体であれば、トリプルでの採取長は80cm程度以上必要でしょう。これは厳しいですね。静的3軸＋動的3軸なら、積算上は2本、実際は3本採ることになるでしょうね（積算にφ75mmの中型3軸はありませんが）。

調査計画を立てるためには、各試験の細かいところまで把握しておく必要がありますね。

単6電池

スタイラスペンの電池を買いに行こうと電気屋さんへ。

AAAAの単6電池なのですが、売っていません。
PC屋かと思い行ってみましたが、これも置いていません。

結局、ネットで購入。まだまだ浸透していないんですね。

3軸CUBと限界状態線

先日、砂質土の三軸CUBについて他の技術者と話をしていたのですが、ちょっと話が合わないところがありました。

p'-q軸上であれば、最悪、1供試体でもφ'を求めることができると考えていたのですが、その方は通常の有効応力経路で数点プロットしないとダメと言われる。

うーん？と思い調べてみると、私の間違いでした。
正確には、正規圧密粘土で c' = 0 とみなせるのであれば、1供試体でも可能です（当然、複数やるべきですが）。しかし、それ以外の土は適用できるかどうか不明。限界状態線の傾き q/p'=M の q/p' を、三軸圧縮試験のモール円のσ'1、σ'3とクーロンの破壊線の傾きφ'で表現し、M = 6sinφ'/(3-sinφ')　の関係を導いてますから。Mを使っていますので、砂は論外ということでしょう。

導出の過程（といっても中学レベルです）は地盤工学会「地盤技術者のためのFEMシリーズ2 弾塑性有限要素法がわかる」 p122です。非排水強度 cu と有効応力のφ'の関係も載っています。が、こちらはモール・クーロンとの違いを理解できていません。拡張ミーゼスを知りませんので。ちょっと調べてみましたが、「拡張」の意図する処がわかりませんでした。

しかし、惜しいですよね。感覚的には砂もCUBであればOKですけど。感覚ですから、ま、通用しません。SYSカムクレイもカムクレイをベースに発展させていますので、いずれは何とかなるかもしれませんが。

Civil3D 2014 とタッチ

出張用PCとして、VAIO Duo 11 SVD1122AJB + Win8 Pro を購入し、1か月がたちました。

http://phreeqc.blogspot.jp/2013/08/vaio-duo-11-svd1122ajb.html

これ、購入して正解でした。

一番良いのは、マウスを使用せずに済むこと。マウスは現場や車内で使えません。PCに搭載されているパッドはマウスに比べ操作性が劣ります（操作が遅くなるのでプチストレス）。そこでペンと指の出番。両手で直感的にタッチ操作できるので、案外手早く操作できます。これだけでも、タッチパネル搭載機種を購入した価値があったと思います。

欠点はスタイラスペンの右クリック。持ち直すたびにボタンの位置を探すか、ペン先で画面を長押しする必要があります。指での長押しは、どこを押しているのか見えないことがあり、細かい作業には向いていません。これらはプチストレス。

あと、現段階で例外的にタッチ操作の難しいのが Civil3D。

2014 からタッチに対応しており、拡大・縮小は指でも可能です。が、マウスほど思うように動きません。また、画面移動が全くできません。2本指でスワイプすれば移動できるはずですが、機能しません。結果的にマウスの使用が必須となりました。惜しい。

あと、「タッチ」パネルも as AutoCAD では表示されるのですが、Civil3D では表示されず。サポートに聞けば不具合とのこと。うーん。

しかし、今回使用してみて、CAD とタッチは相性が良いと感じました。いえ、惜しいと感じることが多かったので、改善に伴い処理速度は格段に向上すると思ったのです。オブジェクト選択後、中を1タップで右クリックメニューが出るようにするとか、ダブルタップにコマンドを追加できるようにしるとか。細かい部分をタップすれば、V-nas のように拡大窓が出るようにしても良いでしょう。まだまだこれからだと思います。

CAD とタッチ、今後が楽しみです。

LSFLOW と LS-RAPID

準三次元地すべり運動解析プログラムによる地すべり性崩壊の被害範囲の予測, 土木研究所資料, 3057, 1992

LSFLOW に使われている支配方程式が掲載されています。それを順に追っていました。

が、最初から？の状態。
いえ、式が難しいわけではなく、掲載されている式に誤りがあるようです。基本的な連続の式の導出から、中の符号、余計なdxがついていたり、必要なところについていなかったり、NがMになっていたり。あちこちに？な部分がありました。ま、こちらも初心者なので、これはこれで意味があるのかもしれませんが。

基本的には、連続の式とニュートンの運動方程式を使っているだけですね。粘性項と外力項があるので、ナビエ・ストークスに近いと思います（自重項はありません）。あとは外力として内部消費が考慮されていますね。
近似はFDM。今まで、なぜ土石流のような大変形が解けるのか不思議で仕方なかったのですが、基本的には変形を解いているのではなく、流体の移動に伴う高さの変化のみを扱っているようです。ですから、粒子が移動するような可視化はできません。納得です。
計算は準三次元とありますが、高さによって物性の変化はないため、2次元に近いですね。

ソースも公開されていますので、それを読んでいますと、ここにも矛盾が。いえ、input ファイルに書かれているパラメーターよりも、読み込むパラメーターのほうが多くなっています。また、地震波の入力や計算のことは記載されていないのですが、ソースは整備されています。作成した時期にズレがあるのでしょうね。ちなみに、ソースは中央開発さん作のようです。
input ファイルの作成は EXCEL ＋ VBA で、結果の可視化は Surfer などで十分対応できそうです。


よく似たソフトに五大開発さんのLS-RAPIDがあります。
作成は佐々先生のようです。論文になっている改良そりモデルでしょうか？
http://www.godai.co.jp/soft/product/products/LS-RAPID/index.htm

こちらも連続の式と運動方程式を使用していますが、粘性をあつかっておらず、自重、外力（せん断抵抗含む）＋αといったシンプルかつ個人的には理解しやすい構成になっています。一応、間隙水圧も考慮されていますが、Jakyの式は3次元応力状態が対象でなかったはずですし、摩擦係数等のフィッティングという点ではLSFLOWと同じ扱いでしょう。離散化も差分、変形でなく高さの変化みの考慮ですので、本質的には同じつくりですね。ま、LSFLOWはすべった後、LS-RAPIDはすべり発生から扱えるという点が、大きな違いかもしれません。
http://phreeqc.blogspot.jp/2011/07/blog-post_30.html

ポスト処理としての動画が掲載されていますが、波動みたいです。うーん。

両者とも、高さのみの変化なので、可視化に工夫が必要なのでしょう。難しいですが。
やはり粒子法でしょうか？こちらはさっぱり進んでいません。理解しないと。

adsorption isotherm と遅延係数

北海道環境保全技術協会が発行している技術レポートを読んでいました。

自然由来重金属や硝酸性窒素といった問題についてのレポートです。数か月前にまとめて購入していたのですが、パッと見、特に目新しいことは書かれていなかったので後回しになっていました。

「吸着工設計マニュアル」に移流分散と、総量の話が出ていました。洩れても基準値を越えなければOKという思想のようです。国もそうですから、ま、良いのでしょう。
http://phreeqc.blogspot.jp/2011/08/blog-post_26.html

現地発生土でしょうか、火山灰土を吸着層へ利用するレポートもありました。この中で adsorption isotherm に関し、Freundlich と Langmuir に対するコメントがありました。
個人的には、ただのフィッティングなんで、線形含めどれでも OK です。ただし、どの isotherm を仮定するかで誘導される式中の R の定義が異なります。
MT3DMS は Visual MODFLOW の移流分散パッケージに採用されており、非線形も含め3種類の isotherm を使用できます。が、Dtransu は線形しか扱えません。例の部長様のように、Dtransu で解くのに Freundlich でフィッティングする、というようななことをしなければ、個人的には何を使っても問題ないと考えています（基礎をおろそかにしてはいけません）。

所詮、遅延係数は計算上のフィッティング係数であり、物理的根拠は薄いと考えます。根拠を求めるなら、反応計算を実施すべきでしょう。
http://phreeqc.blogspot.jp/2010/10/blog-post_31.html

トリプルサンプラー

トリプルサンプラー（ロータリー式三重管サンプラー）をバラした写真です。

硬さが中位以上の粘性土、締まりの程度が中位以上の砂質土を対象としています。

スイベルヘッド。ロッドにつながる側です。ベアリングとバネが入っています。
8個の孔から水が入り、下に抜けます。下側にも同じ数の穴が開いています。
レンチのかかっている管がインナーチューブ。この中にアクリルパイプが入ります。外側にはアウターチューブを被せ、スイベルヘッドにつなげます。水はアウターチューブとインナーチューブの間を流れていきます。インナーチューブの上に1個あいている孔には出口側に逆止弁が入っているため、流入しません。

シューとメタル。先端だけはずしました。86タイプです。

メタルが回っても、ベアリングのおかげでシューは回転しません。

シューがメタルより飛び出していますが、押せばメタルに近づきます（インナーチューブを介し、スイベルヘッド内のバネを押すようになります）。掘削水は孔から出ますが、シューが出ているため試料に接触しません。

試料はシューに押し込まれ、アクリルパイプの中に入っていきます。

これを見れば、砂礫だとシューにあたって掘進できなさそう、、小さな礫であれば入るかも、細粒分がないと落ちそう、などと想像できるのではないでしょうか？

また、シンウォールはシューの部分だけ、デニソンはシンウォール＋メタル＆アウターチューブと考えると理解しやすいですね。孔径を考えると、イマイチ、デニソンの必要性を感じませんが。

先端を外したところ。アクリルパイプ内に試料が収まっています。

いつも思うのが、シューの中の試料を押し込んでよいのかどうか。試験時、先端は使わないので、簡単に押し込めそうなら押し込みます。オペさんに何も言わなければ、100%押し込んでます。これは通常のコア取り時の癖でしょうね。

アクリルパイプを取り出した後、採取長を測り、水抜き、シール、凍結して（細粒分が多いとそのままで）運搬します。

ただ、以前にも書きましたが、今は道具がよくなっていますので通常のφ66ダブルでも乱れの少ない試料は採取可能です。腕も必要ですが。

http://phreeqc.blogspot.jp/2011/06/blog-post_21.html

http://phreeqc.blogspot.jp/2012/10/blog-post_6021.html

その「腕」のナレッジ部分を基準化するのは困難でしょうね。非常に惜しいですが。

ま、お金のかかる取り方が標準かつお客様の御要望であれば、それはそれで仕方ないですが。

ライフジャケット

台船へ行き来するのに、ゴムボートを使用しています。

手漕ぎボートが破れかけたので、エンジン付きのしっかりしたボートを用意していたところ、船上でボーリング作業をしていたオペさんが迎えに来いと合図。

岸にいた一人が迎えに行こうとし、ライフジャケットを付けて手漕ぎボートに乗り込んだ瞬間・・・ドボン。後転しながらひっくり返ってしまいました。一同、騒然。で、大爆笑。

水際の水深50cmだったので爆笑で済んだわけですが、あれが水深数mだと慌てたでしょうねえ。

まず、「携帯大丈夫か！」
次に、「ライフジャケットって、ちゃんと膨らむんだね。」と妙に感心。

今まで、自動膨張のジャケットを着てきましたが、膨らんだところは見たことがありませんでした。落ちたのも初めて見ましたし。あの瞬間を見ると、手で引っ張って膨らませるっていうのは難しいかもしれません。安くて暑いジャケットは浮かないそうですし。

うーん。恐るべしライフジャケット。
重要性を再認識しました。

Dtransu の流速 その2

飽和時の計算について、マニュアルには何も書かれていませんでした。


仕方ないので、ソースを見てみました。

ここですね。


C     DETERMINE RELATIVE HYDRAULIC CONDUCTIVITY
C
        LTB = MTBL3D( MAT )
        DO 61 K=1,NN
          I = LM(K)
C
c-------- confined problem ----

          if( kan2.eq.0 ) then

            CR(I) = 1.0
            TXX(K) = 100.0D0
            go to 61
          end if
c
          IF( P(I).LT.0.0 ) THEN
          CALL INTERP(XYP( 1,NMTB(1,2,LTB)+1 ),TX,P(I),NMTB(2,2,LTB),1)
          TXX(K) = TX
          IF( CR(I).LE.0.0 )
     2    CALL INTERP(XYK( 1,NMTB(1,1,LTB)+1 ),TX,CR(I),NMTB(2,1,LTB),0)
          ELSE
          CR(I) = 1.0
          TXX(K) = 100.0D0
          END IF
C
   61   CONTINUE



      SUBROUTINE SETCR( KODE,P,P1,CR,T,CS,POR,XYK,NK,XYP,NP,XYC,NC )
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      COMMON / CTPRM2/ KAN1,KAN2
C
      DIMENSION XYK(2,NK),XYP(2,NP),XYC(2,NC)
C
      IF( CR.GE.0.0D0 ) RETURN
      IF( KAN2.EQ.0 ) GO TO 2000
      IF( P.LT.0.0D0 )  GO TO 1000
C
C---- SATURATED ---------------------------
C
 2000 CR = 1.0D0
      T  = POR
      CS = 0.0D0
C
      RETURN
C
C---- UNSATURATED -------------------------
C
 1000 IF( KODE.LT.1 ) CM = P
      IF( KODE.GE.1 ) CM = (P+P1)*0.5D0
      CALL INTERP( XYP,T,CM,NP,1 )
      CALL INTERP( XYK,T,CR,NK,0 )
      CALL INTERP( XYC,T,CS,NC,0 )
      T = T*POR*0.01D0
      CS= CS*POR
C
      RETURN
C
C
      END




コメントがないので誤っているかもしれませんが、飽和のフラグを立てると、比透水係数1、体積含水率 = 有効間隙率で計算しているのでしょう。単純に、この条件で負圧となれば不飽和時よりも上昇流や水平方向の流れが発生しやすいということでしょうね。


結局、変換ツールのソースを、圧力水頭がマイナス時は、Vi ⇔ { 0, 0, -0.1} として吐き出すように修正しました。不飽和の部分では鉛直下降流が起こっている、という、よくある簡略化です。

そうすると、TECPLOT の streamline が、思い描いたように綺麗になりました。


うーん、いいのか？

Dtransu の流速

Dtransu の計算結果を Tecplot で可視化しています。

といっても、移流分散ではなく、浸透流。
流速3成分から streamline を描こうとしていたのですが、なかなかできません。

Tecplot の設定は問題なさそうでしたので、変換結果を見ました。すると、節点流速が出ていません。そりゃ、描けないはずです。

Dtransu に付属の変換ツールの詳細を見てみますと、流速を読み飛ばす仕様になっていました。それなら、「読んだ後に節点情報とまとめて書き出せばよい」と考え、ソースを修正。
修正したソースをコンパイルし、実行すると、ちゃんと書き出せました。で、それを Tecplot で読み込み。

が、エラー。
なぜ？と思いながら調べてみますと、原因は私のミスでした。
変換ツールが読み飛ばしていた流速は、要素流速。これを節点流速だと思い込んで、節点情報に付加して書き出していました。Dtransu の浸透流では、要素流速しか書き出さない仕様で、節点流速に直してくれないんですね。移流分散は実流速を節点流速として書き出してくれるのですが。

Tecplot では要素流速の取り込み方がわからなかったのと、それを別途書き出して読み込ませるのも手間なので、結局、ダミーの移流分散として回すことにしました。


計算時間は掛かりましたが、節点流速は書き出されました。
濃度を設定していないため、非数 NaN を 0 に変換し読み込むと、ようやく streamline を描くことができました。長かった。


しかし、これで終わりではありませんでした。
なぜか不飽和の部分で上昇流が発生。今回は飽和の計算しかしていないので、不飽和部分が圧力水頭を持っているとは考えていませんでした。

この続きは、また後日。

斜面災害

8月末から9月頭にかけて、2度の大雨でした。

「大体、2回目の雨で崩れるよね」なんて言ってたら、その日の午後に呼び出し。

国道が崩壊した現場と、市道がすべった現場の2種類でした。出張の合間に、内業を進めようと帰社したら、逆に仕事が増えてしまいました。

国道は完全に地質屋のミス。でも、こういうこともあるのだと、改めて勉強になりました。

市道は早急に観測体制を確立する必要がありました。伸縮計を設置したのですが、通信機や警報機、回転灯などの配線接続は電気工事屋さんみたいな仕事です。もともと、伸縮計の設置は一人ではできなかったのですが、今ではより多くの知識と人数が必要になりました。

当然、その日のうちにネット上で監視できる体制にするため、通信機やサーバーの設定も必要です。今では携帯、スマホで経時グラフを確認できるようになっています。

開発側も大変でしょうが、使う側もそれなりの努力が必要です。

デニソンサンプラーと削孔径

二重管サンプラーの代表として、デニソンがあります。

デニソンサンプラーの中にはシンウォール管が入ります。そのため、掘削径はシンウォールより一回り大きいφ116mmとなります。これは一般的によく知られていますし、正解です。

しかし、二重管サンプラー = φ116mmは誤り。二重管サンプラーもトリプル同様にいくつかの径が出ています。掘削径 86mm で採取できるタイプも販売されています。地盤工学会基準「地盤調査の方法と解説」にもいくつかの径が紹介されています。

しかし、デニソン以外のサンプラーは、あまり見ません。トリプルの86タイプで硬質粘性土も採取できるからでしょうか？

ちなみに、固い粘性土であれば、小径倍圧サンプラーという手もあります。これなら、φ66mmでも採取可能です。ただし、水圧をかけるため、通常のホースでは破裂しますし、場合によっては機械が浮かないようにアンカーを取る必要があります。水周りを耐圧仕様にする必要があるため、やや面倒です。個人的には避けたいサンプラーです。

粘性土の場合、個人的にはアクリル管よりシンウォール管のほうが取り出しやすいので好きです。デニソンで上手に取っていただけるほうがありがたいですね。

トリプルサンプラーと削孔径

ボーリング調査時の乱れの少ない試料採取では、一般的に、以下の削孔径が必要とされています。

シンウォールサンプリング：φ86mm

デニソンサンプリング：φ116mm

トリプルサンプリング：φ116mm

ところがこれ、トリプルのφ116mm は誤りです（積算上はφ116mmしかありませんので正解ですが）。

市販品の中にはφ86mm で採取可能なトリプルサンプラーがあります。昔からよく使われており、普及しています。

私が頻繁に目にするのがコアーパック社の86TRIPLE。

http://www.corepack.co.jp/seihin/packtube_4047/tripletype_4051.html

砂を取るのにφ116mm は必要ありません。φ86mm で十分です。径が小さいと、試料落下のリスクも小さくなります。

オペさんにも何度か話を伺いましたが、サンプラー径が 116mm から 86mm になっても、礫径の影響はほとんど受けないそうです。礫にあたると、どちらも入らない印象が強いようです。それよりも、落下のリスクが小さい、削孔径の小さい方が容易等のメリットを選択される方が多いようです。

メーカーにも話を伺いましたが、一概に「礫が大きければ大きなサンプラーを」、というわけでもなく、「現場に応じて」とのことでした。要は、腕と経験で判断、ということなのでしょう。腕は重要だと思います。

http://phreeqc.blogspot.jp/2011/06/blog-post_21.html

http://phreeqc.blogspot.jp/2012/10/blog-post_6021.html

分業が進み、現場を眺めているだけでは、径の違いにすら気づかないのでしょうね。いまだにトリプル=φ116mm が常識？となっているのはそのせいでしょう。

地震応答の把握

「考え方がよくわかる設計実務 3 耐震設計の基本」大成建設を読んでいます。

これによれば、耐震設計のフローは以下の通りとなっています。

①耐震性能と地震力の設定
②地震応答の把握
③構造物への地震入力と解析
④構造物の健全性照査

②まで読みました。よくわかります。が、これを読むまで②と③が頭の中で混ざっていました。②は地盤がどのように動くか、③は②で求めた地震外力を何らかの形で構造物に作用させ、構造物の変位や断面力がどうなるか、を求める過程ですね。橋梁屋さんや建築屋さんなどにとっては③のほうがメインでしょう。

②については以下の通り。

• 簡易設定法
  設計水平震度、応答変位（地盤の固有周期と速度応答スペクトルを利用）の2種。
• 等価線形解析法
  結果として等価となる G, h を用いる。振動数領域の動的解析法としてSHAKE が代表。
• 逐次非線形解析法
  振動系の運動方程式を時間領域で逐次的に解いていく方法（直接積分法）。各時間ごとに　G, h を変化させ、時々刻々と計算し応答を求める。解法に Newmark のβ法がよく用いられる。
• 有効応力解析
  液状化の検討に使用される。透水を考慮しないFLIP、考慮する（応力と浸透を Biot で連立させて解く） LIQCA が代表。

Newmark のβ法の説明はわかりやすいですね。簡単な積分と図が示されています。なぜβを0.25とすることが多いのかは、今まで「そんなものだ」と思っていたのですが、そうすることで平均加速度（一定）になるんですね。うーん、何でも読んでみるものです。

気象庁の特別警報

今晩から気象庁の特別警報が始まります。

雨に関する発報指標は土壌雨量指数。大雨注意報や警報と同じです。
http://www.jma.go.jp/jma/kishou/know/tokubetsu-keiho/kizyun.html

土壌雨量指数は全国一律のパラメータで計算されています。気象庁で採用された当時、そのパラメータは公開されていませんでしたが（私の調べ方が足りなかったのかもしれませんが）、現在は公開されています。
http://www.jma.go.jp/jma/kishou/know/bosai/dojoshisu.html

計算は全国一律なのですが、得られた結果に対する評価（発報基準値）は地域毎に異なっています。同じ土壌雨量指数が算出されても、大雨警報などで取り扱う場合には、その数値の意味が地域によって異なるということです。詳細は知らないのですが、雨慣れしている地域は基準値を高めに、雨慣れしていない地域は低めに設定されているようです（さすがによくできています）。今回の特別警報では、数十年に1度の値となっていますので、雨慣れというとらえ方は、大外れではないでしょう。論文を調べたら出てくるとおもいますが。
http://www.jma.go.jp/jma/kishou/know/kijun/index_shisu.html

よく似た情報で土砂災害警戒情報があります。これは特別警報や大雨警報の1軸（土壌雨量指数）＋履歴ではなく、多くは2軸（60分積算雨量＋土壌雨量指数）＋履歴のようです。国交省河川局砂防部が詳しい資料を公開されています。
http://www.jma.go.jp/jma/kishou/know/bosai/doshakeikai.html
http://www.mlit.go.jp/river/shishin_guideline/sabo/dsk_tebiki_h1706.pdf



つい先日、先輩の話を聞く機会がありました。
深層崩壊と雨の関係について色々話されていました。目指すところは、土砂災害警戒情報のような2軸指標でなく、１軸の単純な指標だそうです。実際、一連の深層崩壊について実証されていましたので、正解かもしれません。

先輩によれば、素因の話がまだ残っているとのこと。
はい。地質屋は深層崩壊予測に対して明確な答えを出しているレベルにないといえるでしょう。誘因からせめて、説明できない個所の素因を調べる、といった手法のほうが前に進みやすいかもしれませんね。

標高表示板と津波遡上高

最近、街中で標高表示板を見かけるようになりました。


ある街で素人のAさんと、以下のような話をしました。

A「あれは何？」

私「標高表示板ですね。海からどのくらい高いかを示したものです。0.5mですから、海より0.5m高いということです。」

A「津波がきたら、あそこで0.5mの高さになるっていうこと？」

私「いえ、ただ単に、土地の高さを示しているだけです。このあたりだと、津波は○○mまで来ます」

A「それでは全く意味がないじゃないか？あそこでどのくらいの高さまで来るかが示してないと、どこの高さまで逃げればよいかわからないじゃないか？」

絶句。目から鱗でした。
私が表示板を作った訳ではないのですが、なぜか心の中で「すみません」と誤っていました。

確かにそうですよね。私だってよく知らない街への出張で、津波遡上高や避難経路を毎回確認するわけではありません（山は見ますが）。標高表示板の代わりに、津波予想高の表示板のあったほうが、逃げやすいと思います。いくら津波の計算をして、ハザードマップを作っても、伝わらなければ意味がないのです。いえ、そんなことはわかっているつもりでしたが、全くわかっていないことに気づかされました。

やはり、いろいろな方と接しないと、良いものは作れないのかもしれません。

現場透水試験と浸透流解析

「現場透水試験は必要ない。浸透流の逆解析で透水係数を同定すればよいではないか？」

先日、このような相談が私のところに来ました。
極端ですが、ある意味感心しました。

ただ、個人的には反対。
透水係数の妥当性を示す材料（透水試験による透水係数の分布範囲）が欠けると、それはその分布（一連の透水係数を有するとみなせる土質・地質の分布=水理地質）の妥当性を示す材料も同時に欠けるということに繋がります。土質・地質分布の推定誤差を、（分布が正として）透水係数のみに押しつけるようなモデル化は、新たなデータが追加された場合に破綻する可能性が高いのです。バランスが大事でしょう。

初めに設定した透水係数が逆解析結果と大きく異なる原因としては、以下のことが考えられます。
・推定した土質・地質分布が実際と大きく異なる。
・数少ない結果より透水係数が設定されている。（局所的な値が反映されている、粒度のみから推定されている、必要な試験条件を満たしていない、etc）
・同定するデータの種類・量が少ない（計算結果の精度の問題）。

通常は、得られた透水係数の幅を考慮し、モデルをキャリブレーションして行きます。どうしても観測結果を再現できない場合、またはキャリブレーション結果が透水係数の幅や一般値を大きく逸脱する場合、分布の推定を誤っている可能性を疑います。そして、分布を再検討し、計算を続けます。つまり、キャリブレーション作業は、現場で得られた透水係数と推定された土質・地質分布の両者の妥当性を担保する作業だといえるでしょう。

地質屋は現場で観察した根拠をパズルのように組み合わせ、矛盾なく説明できる土質・地質モデルを作成します。しかし、その矛盾のないモデルは大抵複数できてしまいます。その中でも最も可能性の高い（と思う）もの1枚だけを土質・地質モデルとして残します。キャリブレーションがうまくいかない場合は、構築した理屈に沿ってモデルを微修正するか、頭の中に残された2枚目のモデルを試すことになります。それでも合わないなら、計算が合うのはどういった土質・地質分布かを考え、見落とし箇所にあたりを付け、再び現場に戻ります。計算が、自分の推定した土質・地質分布の誤りを指摘してくれるわけです。
これは水の計算だけでなく、変形も同じ事です（水のほうがはるかに単純なので、同視はできませんが）。特に、三次元解析では土質・地質分布が重要ですね。個人的には、通常は透水係数や変形係数、c・φなどの設定誤差より、土質・地質分布の推定誤差のほうがはるかに大きいと感じています。

自然を計算に乗せて解釈したいのであれば、現場をよく歩くこと、原位置試験を実施することなど、基本をおろそかにしないことが近道だと思います。

SYSカムクレイモデル その3

SYSカムクレイの講習会に参加。振り返ると、2年がかりですね。
http://phreeqc.blogspot.jp/2011/11/sys-2.html

講習会参加の許可を貰う際に、上司より、どのように業務に役立つか部長様に説明するよう求められました。
が、面倒なので、適当に流して私費で参加しました。例の部長様に、新しい理論の良さを理解して頂く程の力量も根気もありません。また、「どのように業務に役に立つか？」を知りたがる方は、それが発注の仕様に入らない限り納得されないでしょう。ノルマのあるサラリーマン技術者としては、現時点で流すのがBESTです。

今回は個人の力量を稼ぐべく参加したのですが、私にとってちょうど良い内容でした。
第1回であるためか、カムクレイの限界から、太田・関口モデル、下負荷、そして上負荷、SYSカムクレイと順に話が進みました。以前にFEMや連続体力学の講習等で聞いていたためか、非常に理解しやすい説明でした。

私の解釈も含めて概念を整理すると、以下の通り。（実際に手を動かしたわけではないので解釈に誤りはあると思います。誤りとわかった段階で、今後、修正します。）



①オリジナルカムクレイモデル

• 標準圧密試験（1次元圧密）q=0と等方圧密試験q≠0で正規圧密線（NCL）の傾きは同じ（平行）。しかし、1次元圧密のほうが、体積が小さくなる。p'が同じ（ex. σ’⇔ {1, 1, 1 }と {1.2, 0.9, 0.9} ∴ p’=1） でも、1次元圧密試験では、せん断応力qが生じる。p’によって生じる体積変化が圧密に関する項、qによって生じる体積変化（切片の差）がダイレイタンシー。
• 限界状態線(CSL：Critical State Line)もNCLと平行。
• 練返した土はp’,qが決まれば、vが決まる（排水・非排水せん断で異なるvとならない）。>>>q～p’～v空間で一つの面を形成する。これがロスコー面（Roscoe surface）
• NCLの切片をΝ、CSLをΓとすると、ロスコー面は切片がΝからΓまで連続的に変化する直線の集合ととらえることができる。>>>ΝとΓを線形補間してみよう！・・・これがオリジナルカムクレイモデル。

②-1：修正カムクレイモデル

• オリジナルカムクレイの線形補間を曲線補間にした（ダイレイタンシーのモデル化が異なる、特異点がなくなり使い勝手が良い）。・・・これが修正カムクレイモデル。
• ほかにも、指数関数型（ECモデル）、対数関数型（LCモデル）などあり。

②-2：関口・太田モデル

• オリジナルカムクレイは等方圧密。>>>異方(K0)圧密を考慮してみよう！・・・これが関口・太田モデル（特異点は残るので処理が必要）。
• ほかにも修正カムクレイに異方性を考慮したものもあり。

③下負荷面・上負荷面の追加（自然堆積粘土への適用）

• 上記のモデルは練返し粘土の挙動を説明。自然堆積粘土では以下の挙動が説明できなかった。
• 弾性状態（過圧密）から弾塑性状態（正規圧密）へのなめらかな変遷（過圧密の解消）、繰り返し載荷による体積ひずみの蓄積>>>オリジナルカムクレイや関口・太田モデルにオプションとして下負荷面を追加。
• NCLの外側（限界状態線の外側：impossible state）で大圧縮する（構造の破壊）。>>>さらに上負荷面を追加。修正カムクレイ+異方性（関口・太田、橋口）+下負荷面＋上負荷面=SYSカムクレイ

案外、対処療法的な加筆•修正の繰り返しなのですね。まだまだ発展しそうです。



理解していない点も多くあります。

• 「過圧密の解消」=「土粒子のインターロッキング・ボンドが消失し、塑性膨張が起こる」といった説明は、ダイレイタンシーによる膨張で、構造の高位化を意味する？つまり、正規圧密状態よりもさらに構造の低位化した状態があるということ？
• 砂の締固め（排水）や液状化（非排水）は、構造の破壊（低位化）で説明されています。では、再液状化を再現できるのでしょうか？

また2年後にでも答えを見つけましょう。

Windows 8 とクラウド

VPN は、リモートネットワークのゲートウェイを使わないように、TCP/IPの中のチェックを外せば、 OK でした。以前までは Cisco のクライアントを使用していましたので気が付きませんでした。

ほか、Windows 8 Pro で気になった点。

• メールアプリは POP をサポートしなくなっています。基本、クラウドやウェブベースの考え方なのでしょうね。
• カレンダーアプリが Google カレンダーをサポートしなくなっています(厳密には逆)。gmail calendar といったアプリだと同期可能なようですが、これはもう少し動向を見て、何をベースにするか判断しないといけないでしょう。Google も強くなりすぎていますので。
• アカウントが PC 単体固有（ローカルアカウント）とネット認証（MSアカウント）の2種になっています。後者のほうが利用できる機能が多く、個人的には便利。いままで Win Live や SkyDrive を積極的に利用していたわけではないのですが、今回、live垢をMS垢に切り替えました。

便利になった反面、情報漏洩前提でPCを触っていく気構えが必要かもしれません。

単体からクラウドへ。この変化を Executive 達に浸透させるのは無理でしょうね。

************************************************************
2013.8.17追記
ゲートウェイのチェックを外すと、リモートデスクトップがダメ。うーん。