JP2015141122A

JP2015141122A - 汚染物質の移流拡散解析方法

Info

Publication number: JP2015141122A
Application number: JP2014014549A
Authority: JP
Inventors: 森田　修二; Shuji Morita; 修二森田
Original assignee: Okumura Corp
Current assignee: Okumura Corp
Priority date: 2014-01-29
Filing date: 2014-01-29
Publication date: 2015-08-03
Anticipated expiration: 2034-01-29
Also published as: JP6195377B2

Abstract

【課題】オイラリアン・ラグラジアン手法による移流拡散解析方法において、数値分散が生じるのを抑制して、地下水汚染の精度の良いシミュレーションを行うことができる汚染物質の移流拡散解析方法を提供する。
【解決手段】オイラリアン・ラグラジアン手法による基礎方程式における移流項の解析は、基準要素ごとに所定のタイミングで行われて、各節点の汚染物質の濃度が、移流元の仮想節点の濃度によって更新されるようになっており、基準要素ごとに移流項の解析が行われる所定のタイミングは、1又は２以上の所定のタイムステップの後の、移流元の仮想節点Ａ2から当該基準要素Ｅ0の節点Ａ0までの移動量が、ｘ方向、ｙ方向、又はｚ方向の少なくとも１方向において、メッシュ分割された基準要素の略要素幅に相当する移動量となるタイミングとなっている。
【選択図】図２

Description

本発明は、汚染物質の移流拡散解析方法に関し、特に、移流拡散方程式を解く有限要素解析プログラムを用いた地下水による汚染物質の移流拡散解析方法に関する。

東日本大震災では、自然由来の重金属類を含んだ津波堆積物による土壌汚染が報告されている。今後は、降雨などにより津波堆積物から流出する重金属物質が、地下水汚染を引き起こすことが予想される。また、震災地で発生したガレキの集積場から、人工物の重金属が拡散することも考えられる。このような地下水汚染への対策を考える場合には、広範な領域の地下水汚染のシミュレーションが必要となる。

コンピュータを使用して地下水汚染のシミュレーションを行う方法として、移流拡散方程式を解く有限要素解析プログラムを用いたものが知られている（例えば、特許文献１参照）。また、有限要素法により基本メッシュ間の汚染物質の移流、拡散状況を解析する汚染環境予測方法や、該予測方法を実行させるためのプログラムが開発されている（例えば、特許文献２参照）。特許文献２の汚染環境予測方法では、くぼ地における汚染物質の濃度分布を予測するために、オイラー法として公知の、オイラー型輸送・拡散モデルによる差分法が用いられている。

特開２０１０−６４００２号公報特開２００８−１５１６９４号公報

オイラー法による移流拡散解析方法は、下記の〔数１〕の式に示すような、拡散現象と移流現象とを表す移流拡散方程式を数値解析手法として採用して、タイムステップごとに、〔数１〕の式が解析領域内で釣り合う（濃度が保存される）ように連立方程式を解いて行くものであるが、地下水の流速が早い場合（ほとんどの場合、移流現象を表す流速の方が圧倒的に早い）には、拡散現象と移流現象は、濃度の変化に与える影響度が全く異なり、移流していく中で拡散する汚染物質の濃度は僅かである。このように、汚染物質の濃度の変化に与える影響度が拡散現象と移流現象とで全く異なると、これらを同時に解く手法であるオイラー法では、例えばペクレ数（拡散と移流の比で、拡散係数と、流速×メッシュ幅との比）≒１程度とするために、メッシュ分割を、土粒子の粒径の数倍〜数十倍（数十センチメートル）程度以下のメッシュ幅となるように非常に細かくしなければならないことから、例えば地下水汚染のシミュレーションを広範な領域に亘って行う必要がある場合には、採用することは現実的には不可能である。

一方、このようなオイラー法による移流拡散解析方法の間題点を解決できるようにした解析方法として、オイラリアン・ラグラジアン手法による移流拡散解析方法（ＥＬ法）が知られている。オイラリアン・ラグラジアン手法による移流拡散解析方法は、拡散現象と移流現象とを同時に解くことをせずに、拡散現象と移流現象を順番に解いていくという考え方であり、タイムステップが更新される度に、まず数１の式の第４項〜第６項によって汚染の移動（移流）だけを計算し、次に〔数１〕の式の第１項〜第３項によって汚染の広がり（拡散）を計算するものであり、このように〔数１〕の式を段階的に解くことで、オイラー法の問題点を解消できるようになっている。

しかしながら、有限要素解析プログラムを用いた解析手法では、汚染物質の濃度の値を、メッシュ分割された基準要素の各節点の濃度で表すことになるため、従来のＥＬ法では、タイムステップが更新される際に、移流元の仮想節点（場所によって異なる地下水の流速や流向にしたがって、所定のタイムステップの後に、当該基準要素の各節点まで移動してきた移流元として算定される仮想の節点）が、例えばこれらの仮想節点が含まれる基準要素の中央部分に位置している場合には、移流元の仮想節点の汚染物質の濃度を、例えばアイソパラメトリック要素の形状関数を用いて内挿補間することにより算定して、基準要素の各節点の既知の濃度から求める際に、濃度が急変するような場所では、正しい濃度分布から濃度が分散して、数値分散（数値解析の離散化手法に伴う精度低下）を生じることなる。

このような数値分散を避けるために、従来のＥＬ法では、移流元の仮想節点の位置が、基準要素の各節点の位置に常に近くなるように、クーラン数（流速×タイムステップ間隔と、メッシュ分割された基準要素の要素幅との比）≒１となるように、タイムステップを調整しているが、地下水の流れは、例えば地下水汚染のシミュレーションを行う広範な領域に亘って、一定の流速や流向であることはなく、メッシュ分割された各々の基準要素の領域で、流速や流向が異なっているため、クーラン数が略１となるようにタイムステップを調整したとしても、メッシュ分割されたいずれかの領域では数値分散が生じることで、精度の良いシミュレーションを行うことが困難になる。

本発明は、オイラリアン・ラグラジアン手法による移流拡散解析方法において、数値分散が生じるのを効果的に抑制して、地下水汚染の精度の良いシミュレーションを行うことのできる汚染物質の移流拡散解析方法を提供することを目的とする。

本発明は、移流拡散方程式を解く有限要素解析プログラムを用いた地下水による汚染物質の移流拡散解析方法であって、前記移流拡散方程式として、オイラリアン・ラグラジアン手法による基礎方程式が用いられ、該基礎方程式における移流項の解析は、メッシュ分割された基準要素ごとに、所定のタイミングで行われて、各節点の汚染物質の濃度が、移流元の仮想節点の汚染物質の濃度によって更新されるようになっており、前記移流元の仮想節点は、前記有限要素解析プログラムによる浸透流解析によって解析された、場所によって異なる地下水の流速や流向にしたがって、1又は２以上の所定のタイムステップの後に、当該基準要素まで移動してきた移流元として算定される仮想の節点であり、基準要素ごとに移流項の解析が行われる前記所定のタイミングは、1又は２以上の所定のタイムステップの後の、前記移流元の仮想節点から当該基準要素の各節点までの移動量が、ｘ方向、ｙ方向、又はｚ方向の少なくとも１方向において、メッシュ分割された基準要素の略要素幅に相当する移動量となるタイミングとなっており、基準要素の各節点の濃度が更新される、前記移流元の仮想節点の濃度は、当該移流元の仮想節点が含まれる基準要素の各節点の既知の濃度から、内挿補間することによって算定されるようになっている汚染物質の移流拡散解析方法を提供することにより、上記目的を達成したものである。

そして、本発明の汚染物質の移流拡散解析方法は、基準要素の各節点の濃度が更新される、前記移流元の仮想節点の濃度が、当該移流元の仮想節点が含まれる基準要素の各節点の既知の濃度から、局所座標系のアイソパラメトリック要素の形状関数を用いて内挿補間することによって算定されるようになっていることが好ましい。

また、本発明の汚染物質の移流拡散解析方法は、前記オイラリアン・ラグラジアン手法による基礎方程式における拡散項の解析が、全ての基準要素において、所定のタイムステップごとに行われて更新されるようになっていることが好ましい。

本発明の汚染物質の移流拡散解析方法によれば、オイラリアン・ラグラジアン手法による移流拡散解析方法において、数値分散が生じるのを効果的に抑制して、地下水汚染の精度の良いシミュレーションを行うことができる。

本発明の好ましい一実施形態に係る汚染物質の移流拡散解析方法によって汚染物質の移流拡散が解析される地盤の有限要素モデルを例示する斜視図である。本発明の好ましい一実施形態に係る汚染物質の移流拡散解析方法による解析フローを示すチャートである。浸透流解析による流れ解析の解析フローを示すチャートである。地下水の流れを例示する説明図である。累積移動距離の計算フローを示すチャートである。基準要素と移流元の仮想節点との関係を示す説明図である。局所座標系のアイソパラメトリック要素の形状関数を用いた内挿補間の説明図である。 Warricの実験概要の説明図である。 Warricの実験結果と、従来ＥＬ法及び改良ＥＬ法（本発明の汚染物質の移流拡散解析方法）による解析結果とを比較して示すチャートである。

本発明の好ましい一実施形態に係る汚染物質の移流拡散解析方法は、コンピュータに組み込まれた、移流拡散方程式を解くことが可能な公知の有限要素解析プログラムを用いて作成された、図１に示すような汚染物質の拡散が予想される地盤の有限要素モデル１０を用いて、当該地盤における汚染物質による地下水汚染を、オイラリアン・ラグラジアン手法によって、コンピュータを用いて精度良くシミュレートできるようにする解析方法として採用されたものである。

また、本実施形態の汚染物質の移流拡散解析方法は、従来のオイラリアン・ラグラジアン手法に改良を加えることで、数値分散が生じるのを抑制することにより、〔数２〕に示す公知の（式-１）による浸透流解析と、公知の（式-２）による移流拡散解析とによって、従来のオイラリアン・ラグラジアン手法による解析方法と比較して、汚染物質の拡散が予想される地盤における汚染物質による地下水汚染の状況を、より精度良くシミュレートできるようにする機能を備えている。

そして、本実施形態の汚染物質の移流拡散解析方法は、移流拡散方程式を解く有限要素解析プログラムを用いた地下水による汚染物質の移流拡散解析方法であって、移流拡散方程式として、（式-２）に示すオイラリアン・ラグラジアン手法による基礎方程式が用いられ、（式-２）の基礎方程式における第４項〜第６項の移流項の解析は、例えば図１に示すようにメッシュ分割された基準要素ごとに、所定のタイミングで行われて、各節点の汚染物質の濃度が、移流元の仮想節点の汚染物質の濃度によって更新されるようになっている。

図６に示すように、移流元の仮想節点Ａ2は、有限要素解析プログラムによる浸透流解析によって解析された、場所によって異なる地下水の流速や流向にしたがって、1又は２以上の所定のタイムステップの後に、当該基準要素Ｅ0の節点Ａ0まで移動してきた、上流側の移流元として算定される仮想の節点となっている。

基準要素ごとに移流項の解析が行われる所定のタイミングは、図２のステップＳ3〜Ｓ6に示すように、1又は２以上の所定のタイムステップの後の、移流元の仮想節点Ａ2から当該基準要素Ｅ0の節点Ａ0までの移動量が、ｘ方向、ｙ方向、又はｚ方向の少なくとも１方向において、メッシュ分割された基準要素の略要素幅に相当する移動量となるタイミングとなっている（図６では、Ａ2からＡ1への移動と、Ａ1からＡ0への移動による、２タイムステップ）。

基準要素Ｅ0の節点Ａ0の濃度が更新される、移流元の仮想節点Ａ2の濃度は、当該移流元の仮想節点Ａ2が含まれる基準要素（図６における移流元の仮想節点Ａ2の濃度を算定するための基準要素Ｅ1）の各節点の既知の濃度から、内挿補間することによって算定されるようになっている。

また、同様に、基準要素Ｅ0の節点Ｂ0，Ｃ0，Ｄ0は、それぞれ移流元の仮想節点Ｂ2，Ｃ2，Ｄ2から、仮想節点Ｂ1，Ｃ1，Ｄ1を経て、２タイムステップ後に移動してきた節点とされる。基準要素Ｅ0の各節点Ｂ0，Ｃ0，Ｄ0の濃度が更新される、移流元の仮想節点Ｂ2，Ｃ2，Ｄ2の濃度は、仮想節点Ｂ2については基準要素Ｅ2の各節点の既知の濃度から、仮想節点Ｃ2については基準要素Ｅ0の各節点の既知の濃度から、仮想節点Ｄ2については基準要素Ｅ3の各節点の既知の濃度から、各々内挿補間することによって算定されるようになっている。

本実施形態では、オイラリアン・ラグラジアン手法による基礎方程式を用いて、汚染物質による地下水汚染を解析するのに先立って、図２のステップＳ1の浸透流解析による流れ解析を行うための、例えば図１に示すような、汚染物質の拡散が予想される地盤の有限要素モデル１０を作成する（図３のサブステップＳ1-1〜Ｓ1-4参照）。

汚染物質の拡散が予想される地盤の有限要素モデルを作成するには、図３に示すように、汚染物質の拡散が予想される対象となる地盤の地形データや地層データを入手して、これらのデータをコンピュータに入力する（サブステップＳ1-1）。

すなわち、サブステップＳ1-1では、例えば数値地図データ等から対象となる地盤の地表面高さの情報を入手すると共に、対象となる地盤で実施された土質調査や揚水試験の結果から、深度ごとの土質や透水係数等の、地層ごとの透水性に関する情報を入手して、これらのデータをコンピュータに入力する。これらのデータは、数値解析を行うために、メッシュ分割した６面体等の要素で表現された、対象となる地盤の有限要素モデル１０を作成するための基本的なデータとなる。

次に、サブステップＳ1-2で、現地における地下水位に関する情報を入手する。すなわち、サブステップＳ1-2では、例えば気象庁や海上保安庁、或いは国土交通省等から公開されている、対象となる地盤を含む現地の降雨量や海岸線の潮位変動、河川水位等の情報を入手すると共に、数値解析によって対象となる地盤中の地下水の流れを再現するための、地下水を涵養する水源のデータや土質調査による地下水位のデータ等を入手して、これらのデータをコンピュータに入力する。

次に、サブステップＳ1-3で、汚染物質の拡散が予想される地盤の解析領域を設定する。すなわち、サブステップＳ1-3では、対象となる地盤を含む現地において、汚染源１１（図１参照）を中心に汚染が広がる可能性のある地域を想定すると共に、地下水が上流から下流へと汚染が広がることを考慮して、地下水流の影響範囲や対象となる汚染源１１による汚染範囲から、解析領域を設定する。解析領域は、例えば汚染源１１を中心とする数百メートル四方の領域として設定することが可能である。

そして、サブステップＳ1-4では、対象となる地盤を有限要素モデル１０としてモデル化した、メッシュデータを作成する。すなわち、サブステップＳ1-4では、例えばサブステップＳ1-1及びサブステップＳ1-3においてコンピュータに入力された情報に基づいて、コンピュータに組み込まれた公知の有限要素解析プログラムの機能によって、対象となる地盤をモデル化したメッシュデータによる有限要素モデル１０を作成する。メッシュ分割（要素の分割）にあたっては、汚染源１１の周辺や、地下水の流速が大きい場所、地下水の流況が変化しやすい場所等では、特にメッシュ分割を細かくして、汚染濃度の急激な変化や地下水の流況を、精度よく表現できるようにすることが可能である。

本実施形態では、上述のサブステップＳ1-1〜Ｓ1-4によって、例えば５００ｍ×４００ｍの領域を解析領域として、節点数が４８×３９×２１＝３９３１２、要素数が４７×３８×２０＝３５７２０、地表面付近の鉛直分割をΔｚ＝０．４ｍとした、例えば図１に示すような有限要素モデル１０が形成さる（図１では、汚染源１１の周囲の部分が抽出されて描かれている。）。

また、本実施形態では、サブステップＳ1-5で、タイムステップの設定が行われる。有限要素解析プログラムによる数値解析では、空間的にはメッシュ分割した格子点や要素によって事象を表現し、時間的にはタイムステップで設定した時刻ごとに繰り返し計算を行って事象を表現することになる。サブステップＳ1-5では、例えば従来のＥＬ法と同様に、地下水の流速の変化等を考慮すると共に、サブステップＳ1-4で作成したメッシュ分割に留意して、可能な限り数値分散が抑えられるように、タイムステップを調整して設定する。サブステップＳ1-5で設定されたタイムステップは、全ての解析領域において、同一のタイムステップが適用される。サブステップＳ1-5では、事象の変動を考慮して、変化が激しい段階では細かくタイムステップを分割することが好ましい。タイムステップは、例えば数分〜数日の範囲で設定されることが好ましい。例えば、図１に示す有限要素モデル１０において、降雨量を年間の平均の降雨量である約５ｍｍ／日、海岸線の潮位変動を１日周期と設定した場合には、１タイムステップを０．２５日とすることができる。豪雨時などに降雨量が１０倍となった場合には、１タイムステップを０．２５日の１／１０とすることが好ましい。タイムステップは、各節点の地下水位（水圧）が大きく変動しないように、例えばメッシュ分割や水源データを考慮して、適宜設定することができる。浸透流解析と移流拡散解析とは、同一のタイムステップで逐次計算を行うのが通常であり、多くの場合は浸透流解析の方が適切なタイムステップが細かくなるため、移流拡散解析での数値分散を助長することになるが、本実施形態では、数値分散を意識せずに、浸透流解析に適切なタイムステップを設定することが可能になる。

さらに、本実施形態では、サブステップＳ1-6で、設定されたタイムステップごとに、上述の〔数２〕の（式-１）の浸透流方程式を解くことで、地下水の流れを算定する。すなわち、サブステップＳ1-6では、上述のサブステップＳ1-1〜Ｓ1-5によって得られたデータに基づいて、有限要素解析プログラムによって、（式-１）の浸透流方程式を解くことで、地盤をモデル化した有限要素モデル１０のメッシュデータにおける、各節点の地下水位（水圧）を求めることができ、これによって地下水の流れを算定することが可能になる。

さらに、本実施形態では、サブステップＳ1-7で、設定されたタイムステップごとに、有限要素モデル１０のメッシュデータの全てについて、要素ごとの流速と流向を算定する。すなわち、サブステップＳ1-7では、サブステップＳ1-6で求めた各節点の地下水位（水圧）に基づいて、有限要素モデル１０のメッシュデータの全要素について、流速と流向を算定する。算定された流速と流向を用いることによって、後述する上述の〔数２〕の（式-２）による移流拡散解析において、流れによる汚染物質の移流が計算される。地下水の流れは、地形や地層、水源によって影響を受けるので、有限要素モデル１０のメッシュデータの要素ごとの流速や流向は様々で、日々刻々と変化することは明らかであり、従来のＥＬ法では数値分散が大きくなる原因となっていた。本実施形態では、後述するように、（式-２）による移流拡散方程式の移流項の解析を、所定の条件下で行うことで、数値分散が大きくなるのを効果的に抑制できるようになっている。

図４は、流速と流向が変化する例を示したもので、矢印の長さと向きが流速と流向を表している。下段の要素ほど流速が低下しており、この場合、要素の大きさに対して相対的に流速が小さいので、従来のＥＬ法による移流拡散解析では、下段の要素ほど数値分散が大きくなると考えられる。

本実施形態では、ステップＳ1で浸透流解析による流れ解析を行なったら、サブステップＳ1-4で作成された有限要素モデル１０のメッシュデータを用いて、図２に示すように、ステップＳ2で、（式-２）の移流拡散方程式（基礎方程式）による移流拡散解析における拡散項の解析を行う。すなわち、ステップＳ2では、従来のＥＬ法と同様に、（式-２）の移流拡散方程式における第１項〜第３項の拡散項の解析が、コンピュータに組み込まれた有限要素解析プログラムによって、全ての要素（基準要素）において、サブステップＳ1-5で設定された所定のタイムステップごとに行われて更新されるようになっている。

そして、本実施形態では、図２のステップＳ3〜Ｓ6に示すように、（式-２）の移流拡散方程式における、第４項〜第６項の移流項の解析が行われる。図２のステップＳ3〜Ｓ6の移流項の解析では、全要素について所定のタイムステップごとに行われていた従来のＥＬ法による移流項の解析とは異なり、有限要素モデル１０を構成する基準要素の各々について、基準要素ごとに移流項の解析が行われるようになっており、所定のタイムステップごとに、後述する所定の条件を満たした基準要素についてのみ解析が行われて、解析が行われた基準要素の各節点についてのみ、汚染物質の濃度が移流元の仮想節点の濃度によって更新されるようになっている。

すなわち、図２のステップＳ3〜ステップＳ6による移流項の解析では、所定のタイムステップごとに、サブステップＳ1-6やサブステップＳ1-7によって得られたデータに基づいて、移流による移動距離の計算がなされ（ステップＳ3）、移動距離が基準要素の基準幅よりも大きいか短いかを判定する（ステップＳ4）。ｘ方向、ｙ方向、及びｚ方向の全てにおいて移動距離が基準幅よりも短い場合には、移流項の解析を行うことなく、移動距離だけを記憶して、次のタイムステップに移り、ｘ方向、ｙ方向、又はｚ方向のいずれかににおいて移動距離が基準幅よりも大きくなっている場合には、移流項の解析を行って（ステップＳ5）、移動距離をリセットするようになっている（ステップＳ6）。

ここで、移流による移動距離の計算するステップＳ3では、図５のサブステップＳ3-1及びＳ3-2に示すように、ｘ方向、ｙ方向、又はｚ方向における、移流による累積移動距離を計算するようになっている。すなわち、サブステップＳ3-1では、１タイムステップごとに、ｘ方向、ｙ方向、及びｚ方向への移動距離の計算がなされると共に、サブステップＳ3-2では、複数のタイムステップによる、ｘ方向、ｙ方向、及びｚ方向への累積移動距離の計算がなされるようになっている。

図５において、Δｔは１タイムステップでの時間の長さであり、１タイムステップで汚染物質が移動する距離は、流速と時間の長さの積（Ｖ×Δｔ）で表されることになる。流向は流速のｘ、ｙ、ｚ成分（Ｖ_X 、Ｖｙ、Ｖｚ）で表されるので、移動方向は移動距離のｘ、ｙ、ｚ成分（ΔＬ_X、ΔＬy、ΔＬz）で表されることになる。連続する複数のタイムステップの間に移動する累積移動距離（Ｌ_X、Ｌy、Ｌz）は、タイムステップ毎の移動距離の和（ΣΔＬ_X、ΣΔＬy 、ΣΔＬz ）で表されることになる。

移動距離又は累積移動距離が基準要素の基準幅よりも大きいか短いかを判定するステップＳ4では、移動距離又は累積移動距離と基準要素の大きさの相対関係をチェックする。すなわち、ｘ方向、ｙ方向、及びｚ方向における移動距離又は累積移動距離（Ｌ_X、Ｌy、Ｌz）が、基準幅よりも大きくなっているか否かを判定する。そして、移動距離又は累積移動距離（Ｌ_X、Ｌy、Ｌz）のいずれかが基準幅（Ｒ×ｄx、Ｒ×ｄy、Ｒ×ｄz）より大きくなった場合（Ｌ_X＞Ｒ×ｄx、Ｌy＞Ｒ×ｄy、Ｌz＞Ｒ×ｄz）に、移流項を計算するステップＳ5に進むことになる。Ｒは、基準幅を決定する際に入力される、ほぼ１に近い値となっている係数であり、１よりも小さい値であることが好ましい。

図６は、基準要素と移流元の仮想節点との関係を２次元座標で例示する説明図であり、移流元の仮想節点Ａ2は、有限要素解析プログラムによる浸透流解析によって解析された、場所によって異なる地下水の流速や流向にしたがって、1又は２以上の所定のタイムステップ（本実施形態では２タイムステップ）の後に、当該基準要素Ｅ0の節点Ａ0まで移動してきたと考えられる、上流側の移流元として算定される仮想の節点となっている。図６の例では、１タイムステップでは、基準要素Ｅ0の節点Ａ0には、移流元の仮想節点Ａ1から移動してきたことになるが、移流元の仮想節点Ａ1から基準要素Ｅ0の節点Ａ0までの移動距離（ΣＶ_X×ｄt 、ΣＶｙ×ｄt）が、ｘ方向及びｙ方向のいずれにおいても、基準幅（Ｒ×ｄx、Ｒ×ｄy）に満たないため、移流項の計算はまだ行わない。

一方、２タイムステップでは、基準要素Ｅ0の節点Ａ0には、移流元の仮想節点Ａ2が、１ステップ前の移流元の仮想節点Ａ1を介して移動してきたことになるが、移流元の仮想節点Ａ2から基準要素Ｅ0の節点Ａ0までの累積移動距離（Ｖ_X×ｄt 、Ｖｙ×ｄt）が、ｘ方向やｙ方向において、基準幅（Ｒ×ｄx、Ｒ×ｄy）を超えることになるので、移流元の仮想節点Ａ2が、基準要素Ｅ0の節点Ａ0まで移動する２タイムステップのタイミングで、移流項を解析するステップＳ5において、（式-２）の移流拡散方程式における第４項〜第６項の移流項の計算が、コンピュータに組み込まれた有限要素解析プログラムによって行われる。

本実施形態では、上述のように、所定のタイムステップごとに、全ての基準要素について移流項の解析を行うのではなく、移流項を解析するステップＳ5において、移流元の仮想節点から基準要素の各節点への移動距離（累積移動距離）が、ｘ方向、ｙ方向、及びｚ方向のいずれかにおいて、基準幅（Ｒ×ｄx、Ｒ×ｄy、Ｒ×ｄz）を超えたか否かの判定が、基準要素ごとに行われて、基準幅（Ｒ×ｄx、Ｒ×ｄy、Ｒ×ｄz）を超えた基準要素についてのみ、所定のタイムステップごとに、移流項の解析が行われるようになっている。

また、本実施形態では、移流項の解析を行う際に、移流元の仮想節点は、基準要素の各節点における、更新されるべき汚染物質の濃度を算定するために利用される。図６の例では、基準要素Ｅ0の節点Ａ0の濃度は、移流元の仮想節点Ａ2の濃度が移動してきたものと考えられることから、移流元の仮想節点Ａ2の濃度が、更新されるべき基準要素Ｅ0の節点Ａ0の濃度として算定される。移流元の仮想節点Ａ2の濃度は、当該移流元の仮想節点の濃度を算定するための、仮想節点Ａ2が要素内に含まれる、各節点の汚染物質の濃度が既知の例えば隣接する濃度算定用の基準要素Ｅ１から、内挿補間することによって算定することができる。本実施形態では、移流元の仮想節点の濃度は、当該移流元の仮想節点が含まれる基準要素Ｅ１の各節点の既知の濃度から、局所座標系のアイソパラメトリック要素の形状関数を用いて内挿補間することによって、算定されるようになっている。

すなわち、アイソパラメトリック要素の形状関数は、図７に示すように、全体座標系における基準要素Ｅ1の内部の点Ｃ0の濃度を、Ｃ1〜Ｃ4の濃度により図７の右欄の計算式によって表すものである。図７の右欄の計算式のＮ1〜Ｎ4が、アイソパラメトリック要素の形状関数である。全体座標系の基準要素Ｅ1を、局所座標系の基準要素Ｅ1’に写像することで、Ｃ0を内挿するための座標は、局所座標系のr、tに変換することが可能である（マトリックス有限要素法、O.C.ツエンキーヴィッツ著、吉織雅夫・山田嘉昭監訳、培風館参照）。

本実施形態では、ステップＳ4で移動距離が基準要素の基準幅よりも大きいか短いかが判定されて、ステップＳ5において、ｘ方向、ｙ方向、又はｚ方向のいずれかにおいて移動距離が基準幅よりも大きくなった基準要素について、移流項の解析を行ったら、当該移流項の解析を行った基準要素の移動距離が、ステップＳ6においてリセットされて、次のステップに移行すると共に、タイムステップが更新されることになる。

そして、上述の構成を備える本実施形態の汚染物質の移流拡散解析方法によれば、オイラリアン・ラグラジアン手法による移流拡散解析方法において、数値分散が生じるのを効果的に抑制して、地下水汚染の精度の良いシミュレーションを行うことが可能になる。

すなわち、本実施形態では、オイラリアン・ラグラジアン手法による基礎方程式における移流項の解析は、基準要素ごとに所定のタイミングで行われて、各節点の汚染物質の濃度が、移流元の仮想節点の濃度によって更新されるようになっており、基準要素ごとに移流項の解析が行われる所定のタイミングは、1又は２以上の所定のタイムステップの後の、移流元の仮想節点Ａ2から当該基準要素Ｅ0の節点Ａ0までの移動量が、ｘ方向、ｙ方向、又はｚ方向の少なくとも１方向において、メッシュ分割された基準要素の略要素幅に相当する移動量となるタイミングとなっている。

これによって、本実施形態では、移流による移動距離（累積移動距離）と要素幅とを、基準要素ごとに比較判定して、一定の移流がある段階で基準要素ごとに移流項の計算を行うので、場所による流速や流向の相違を移流項の計算に考慮することが可能になると共に、移流元の仮想節点の位置が、基準要素の各節点の位置に常に近くなった状態で、移流項を計算することが可能になるので、数値分散が生じるのを効果的に抑制して、地下水汚染のシミュレーションを精度の良く行うことが可能になる。

図８は、トレーサ水の浸透について、Warricが実施した実験の概要を示すものであり、図９は、Warricが実施した実験の実験結果と、従来ＥＬ法及び改良ＥＬ法（本発明の汚染物質の移流拡散解析方法）による解析結果とを比較して示すチャートである。Warricが実施した実験は、トレーサ水を２．８時間かけて浸潤させた後に、原水を１７．５時間かけて浸潤させ、地表面からの浸潤水の濃度を測定して、鉛直方向の濃度分布の時間的変化を示したものとなっている。実験では６ｍ×６ｍの平面領域にトレーサ水を浸潤させているが、従来ＥＬ法や改良ＥＬ法による解析では、６ｍ×６ｍの中心線上の鉛直断面における２次元問題として簡略化した、平面問題として取り扱っている。

Warricによる実験結果と、従来ＥＬ法や改良ＥＬ法による解析結果とを比較して示す図９によれば、従来ＥＬ法では数値分散による濃度低下が顕著に見られるが、改良ＥＬ法では、数値分散が抑制されていて、Warricによる実験結果のトレーサ水の濃度の浸潤深さが、精度良く評価できていることが判明する。

なお、本発明は上記実施形態に限定されることなく種々の変更が可能である。例えば、基準要素の各節点の濃度が更新される、移流元の仮想節点の濃度は、局所座標系のアイソパラメトリック要素の形状関数を用いて内挿補間することによって算定する必要は必ずしもなく、公知の有限要素解析プログラムによって解析可能な、その他の種々の関数を用いて算定することもできる。

１０汚染物質の拡散が予想される地盤の有限要素モデル
１１汚染源

Claims

移流拡散方程式を解く有限要素解析プログラムを用いた地下水による汚染物質の移流拡散解析方法であって、
前記移流拡散方程式として、オイラリアン・ラグラジアン手法による基礎方程式が用いられ、
該基礎方程式における移流項の解析は、メッシュ分割された基準要素ごとに、所定のタイミングで行われて、各節点の汚染物質の濃度が、移流元の仮想節点の汚染物質の濃度によって更新されるようになっており、
前記移流元の仮想節点は、前記有限要素解析プログラムによる浸透流解析によって解析された、場所によって異なる地下水の流速や流向にしたがって、1又は２以上の所定のタイムステップの後に、当該基準要素の各節点まで移動してきた移流元として算定される仮想の節点であり、
基準要素ごとに移流項の解析が行われる前記所定のタイミングは、1又は２以上の所定のタイムステップの後の、前記移流元の仮想節点から当該基準要素の各節点までの移動量が、ｘ方向、ｙ方向、又はｚ方向の少なくとも１方向において、メッシュ分割された基準要素の略要素幅に相当する移動量となるタイミングとなっており、
基準要素の各節点の濃度が更新される、前記移流元の仮想節点の濃度は、当該移流元の仮想節点が含まれる基準要素の各節点の既知の濃度から、内挿補間することによって算定されるようになっている汚染物質の移流拡散解析方法。
基準要素の各節点の濃度が更新される、前記移流元の仮想節点の濃度は、当該移流元の仮想節点が含まれる基準要素の各節点の既知の濃度から、局所座標系のアイソパラメトリック要素の形状関数を用いて内挿補間することによって算定されるようになっている請求項１に記載の汚染物質の移流拡散解析方法。
前記オイラリアン・ラグラジアン手法による基礎方程式における拡散項の解析は、全ての基準要素において、所定のタイムステップごとに行われて更新されるようになっている請求項１又は２記載の汚染物質の移流拡散解析方法。