パッシブ露凝縮器の水分凝縮

はじめに熱移動モデルの設定
パッシブ露凝縮器モデルの概要凝縮モデル
物理分野の結合付録
気象パラメータ用語集
流体力学モデルの設定参考文献

はじめに

パッシブ露凝縮器(PDC)とは電源を使わずに大気から水分を集める装置である[Beysens, 2018]．この装置は別の水源を提供することを目的としており，科学分野のコミュニティでは，水分を集めることにおけるPDC装置の効率を最適化する方法を探っている．

ここで作成するPDEモデルはモロッコのミルレフトに実際に設置されている装置[Lekouch et al., 2012]と比較し，約24時間でどれほどの水が溜まったかを示す．

この装置の可視化と実験的設定は[Lekouch et al., 2012]で見ることができる．このシミュレーションでは[Flores-Castillo, 2022]のAppendix 1で与えられている，実験で得られたデータを使い，測定された水の量とシミュレーションの結果を比較する．

モロッコのミルレフトを可視化する：

このチュートリアルでは，露水の採取量を推定するために，実際の気象データを使ってPDCの過渡モデルのシミュレーションを行う方法を説明する．このモデルは，材料や傾斜角度の他，周囲温度，相対湿度，雲量，太陽放射等の気象パラメータを含む外部変数等のPDCの表面特性を考慮する．このモデルはよりよいPDC装置を最適化したり設計したりするためのツールとして使うことができる．ここで提示する科学的モデルは，[Flores-Castillo, 2022]に基づいた，PDCの二次元モデルである．

パッシブ露収集器は，露という気象現象のもとで有効である．露の形成は，大気中の水蒸気が表面に凝縮する過程[Beysens, 2018]に対応する．露が表面に形成されるためには，表面の温度は露点より低くなければならない．露点は湿った空気が飽和するように冷えなければならない温度として定義される．つまり，一定圧力と仮定すると，露点では水蒸気は蒸発するのと同じ速度で液体に凝縮するのである．したがって，露点より低い温度の表面では，水は蒸発するより凝縮する量の方が多い．この過剰な水が収集できる．

収集器が熱エネルギーを自然に解放することによって，表面は露点より冷えることができる．また，気象効果と技術的条件もPDC表面の冷却に影響する．

ここで使用するPDC装置は水平面から30°度の角度が付いた平面PDCである．PDCは主に断熱材料の上部の凝縮器表面でできている．凝縮器表面の特性（放射率や熱容量等）は，どれだけの水が凝縮するかに影響するので，非常に重要である[Flores-Castillo, 2022]．

このモデルでは，さまざまな物理現象を考慮に入れる．例えば，周囲空気の速度，空気およびPDC装置自体における熱移動，さまざまな気象パラメータ等である．このすべてが時間とともに変化する現象である．このモデルは複雑な性質のため，実行時間とメモリの両方において，計算的に高価なものになる．

有限要素パッケージをロードする：

パッシブ露凝縮器モデルの概要

このモデルでは，流体力学モデルと熱移動モデルが連続して実行される．その結果が凝縮器モデルに送られる．次のセクションではサブモデルの概要を説明する．

まず，凝縮器と断熱材からなる，パッシブ露収集器とその構成要素を示したスケッチを見てみる．空気がPDCを取り巻いている．空気は左境界から入り，右境界から出ていく．ここでは - 座標系の2D領域を使う．

5.gif

パッシブ露収集器とその構成要素．周囲は空気．

流体力学モデル

流体力学モデルは，パッシブ露収集器の上を流れる空気のシミュレーションを行う．使用する流体力学モデルは，時間依存のナビエストークス方程式である．ベクトル値のナビエストークス方程式は以下で与えられる：

ここで は密度， は粘度， は圧力， はベクトル値の速度場である．速度場のコンポーネントは であり，それぞれ 方向と 方向である．解く従属変数は速度のコンポーネントと相対圧力である．

流体力学モデルの結果は，装置の上を流れる空気の速度場 である．この 流体速度場の解は，PDC装置を取り巻く空気の温度分布をもたらすため，凝縮する水分量を推定したい場合に知る必要がある．

熱移動モデル

PDCシミュレーションの2つ目のコンポーネントとして，熱移動モデルが凝縮器表面の温度を得るために使われる．表面温度を知ることで，凝縮器表面とその周囲の空気の間の熱移動係数 が計算できる．これは表面で凝縮する水の量を計算するのに必要である．

熱方程式は永代領域とPDC装置の両方でアクティブである．過渡熱移動モデルでは，温度分布 は過渡熱方程式で記述される．過渡熱方程式は以下で与えられる：

ここで は定圧熱容量， は温度， は熱伝導率である．上で説明した，流体力学モデルを使って計算した空気流は，装置上で対流する熱をモデル化するのに使われる．それぞれ 方向， 方向のコンポーネント ， を持つ流体流モデルで計算した空気流速度 は対流項 を活発にする．

PDCはさまざまな効果を介して周囲と熱を交換する[Flores-Castillo, 2022]．この装置は対流 によって熱エネルギーを失い，放射 によってエネルギーを失うか吸収するかし，太陽放射照度 および凝縮 によってエネルギーを吸収する．

40.gif

PDCの凝縮器表面における熱移動のさまざまなプロセスの概要．対流項 は装置上の空気流速度をモデル化し，PDCからの熱エネルギーを排除する．放射 は放射によるPCDと周囲との熱エネルギーの流入，流出をモデル化する．流入熱流束と考えられる凝縮プロセス も，装置の温度に影響を及ぼす．流入熱流束と考えられる太陽放射照度 も装置の温度に影響を与える．

凝縮モデル

流体流モデルと熱移動モデルは凝縮モデルに付随して起こる．凝縮モデルによって，凝縮の量が計算できる．凝縮はPDEとしてモデル化されないが積分でモデル化される．

凝縮モデルは，凝縮水の質量変化率 を表す単独の方程式を基にしており，次の方程式で記述される[Beysens, 2018]：

ここで はさらされている凝縮器表面の表面積， は凝縮器の上の湿った空気の水圧を周囲温度 の関数として表したもの， は凝縮器温度 における水蒸気の圧力である．は水蒸気移動係数であり熱移動係数 に比例する．これは以下の方程式で与えられる：

ここで は，水が液体から気体に位相を変化させるのに必要な，蒸発の潜熱である． は空気中の水蒸気の分圧と空気温度を関連付ける湿度定数であり，気象学や気候学等の分野で広く使われている．

方程式(1)を解くためには，これまでに述べた変数のすべて，特に熱移動係数 を知る必要がある．

物理分野の結合

原理上は，凝縮水の量は連続アプローチと分離アプローチの2つの方法で計算できる．この2つの方法の一番の違いは，連続アプローチではPDCの温度は凝縮項 に依存していないということである．凝縮項を考慮しないことによって，モデルを連続して実行できるようになる．つまり凝縮モデルと熱移動モデルは完全には結合していないということである．このため連続アプローチの方が簡単で計算も早くできるが，分離アプローチの方が正確である．

連続アプローチ

連続アプローチではナビエストークス方程式(2)を先に解いてから熱方程式(3)を解く．その後で熱移動係数 と方程式(4)の残りの変数を計算して，水量を算出する．

この連続解法の手順を以下の図に示す．

64.gif

水凝縮を解く連続アプローチ．

気象パラメータ

PDEモデルの入力として実際の気象データを使うことは解析にとってプラスになる．実際の気象データは特定の場所におけるモデルからより現実的で正確な結果を得る役に立つ．さらに実際のデータを使うことで，通常または極端な状況に対する装置の設計を最適化するのにも役立つ．

シミュレーションで使う気象データには周囲温度 ，相対湿度 ，露点 ，風速 ，雲量率がある．雲量率は特定の場所で空が雲に隠れている平均の割合を言う．気象パラメータはすべて時間依存である．

ここで提示するモデルは論文[Flores-Castillo, 2022]のAppendix Iで示されているデータを使う．これらのデータ点は連続関数に変換する必要がある．そのために，データを補完する．

さらに付録セクションでは，WeatherDataで利用可能な世界中の標準的な測候所からの現在および過去の気象データを使う方法を説明する．これを使うことで，世界のどの地点におけるパッシブ露収集器の性能でも予測可能になる．

取得した気象変数はシミュレーション全体で使う．

気象関数

まず使用する気象関数を定義する．この関数には露点 ，天空温度 および ，天空放射率 および ，雲量率が含まれる．

露点 は周囲温度 および相対湿度に比例し，以下で与えられる：

ここでは以下で表せる：

方程式(5)では，周囲温度 は単位で与えなければならない．

方程式5に従って露点 を定義する：

露点 （単位は）および高度 （単位は）に依存する上空放射率 は，以下で与えられる：

方程式7に従って上空放射率 を定義する：

雲量を含む上空放射率は以下で与えられる：

以下の方程式によると，上空温度 は周囲温度 に関係している：

したがって，雲量 を含む上空温度は上空放射率 を に置き換えることで導くことができる：

方程式9に従って上空温度 を定義する：
方程式8および10に従って，雲量の影響のある放射率 と上空温度 を定義する：

後で気象データを使って，周囲温度 ，相対湿度，風速，雲量率の関数を得る．

モロッコのミルレフトからの気象データ

PDCシミュレーションを実行するために，モロッコのミルレフトからの気象データ[Lekouch et al., 2012]を使う．これらのデータは2007年の9月日から9月日までのものである．この場所からの実験データが存在することが分かったら，モデルを比較することができる．

地理パラメータは緯度，経度である．

緯度と経度を定義する：
測地位置を定義する：
高度をで得てに変換する：
開始日をDateObjectとして定義する：
終了日をDateObjectとして定義する：
両方の日付の間の経過時間をで得る：
GeoGraphicsで場所を可視化する：

次にデータをインポートする．

データを抽出するファイル名と行のリストを定義する：
ImportおよびInterpolationを使って，ファイルからデータをインポートし，それを補間関数にする：

風速と雲量には，データを平滑化するためにGaussianFilterが適用された．

残念ながら実際のデータには，負の雲量率等多くの問題がある．このような欠点を解決するために，負の値を処理するためのクリーンアップ関数を作成する．

UnitStepを使って負の値を修正する：

風速データでは，最低風速を指定する必要がある[Flores-Castillo, 2022]．この場合はである．

if文を使って，関数も作成する：

一般に風速は地上の高さで観測される．それより低い高さでの速度が知りたい場合は補外法が必要である．これには通常，以下のような従来の対数変動が使われる．

ここで は風速が理論上ゼロになる高さであり，は地上 の地点での風速である． は一般的に0.1となる[Beysens, 2018]．

したがって，地上のところの風速が知りたければ，方程式(6)に を適用するだけである．

地上 での速度を定義する：
データの長さを得る：
周囲温度 ，天空温度 ，露点 ，雲の影響を受けた天空温度 をプロットする：
雲の影響がある上空放射率 と影響のない上空放射率 の差分を可視化する：
相対湿度と雲量率を可視化する：
地上とにおける風速をプロットする：

流体力学モデルの設定

はじめに

このセクションでは流体力学モデルを設定して解く．PDC装置を取り巻く過渡流れ場 について解くためにナビエストークス方程式を使う．

ベクトル値のナビエストークス方程式は以下で与えられる：

ここで は密度， は粘度， は圧力， はベクトル値の速度場である．その成分は 方向 方向それぞれで である．解く従属変数は速度の成分と相対圧力である．

簡単にするために重力等の外部の力はないものとする．このシミュレーションでは，空気は非圧縮，つまり密度は一定であるとする．また，数値解の複雑性を軽減するために自然対流は無視し，強制対流だけを考慮する[Flores-Castillo, 2022]．

ナビエストークス方程式は，形状の空気領域でのみアクティブである．

流体力学PDEモデルの他，境界条件を指定する必要がある．

  • 流入速度プロファイルをPDC装置の左側に設定する．
  • PCDの壁面と地表レベルに滑りなし境界条件を設定する．滑りなしとは，表面においてすべての方向の速度はゼロということである．これは壁面をモデル化する際に流体力学において一般的な境界条件であるが，PDC表面で滑りなし条件を選ぶことは，以下に示す熱移動モデルでは必須であり，それについては後に説明する．
  • 圧力 が任意の値0に設定されている，領域の右側における流出条件 ．

相対圧力は，基準圧力に関する流体の圧力である．ほとんどの場合，基準圧力 は周囲圧力に設定される．必要に応じて，以下の方程式を使ってシステムの絶対圧力 を計算することができる：

この方程式は，流出条件が1気圧の絶対圧力に設定されていることを示す．初期条件と境界条件は相対圧力値を使って定義しなければならない．

領域

以下の図に示すように，解析の領域は灰色で描いたPDC装置と青で描いた周囲の流体領域である．

160.gif

PDC装置を灰色で，周囲の空気を青で描いた解析の領域Ω．

PDC装置周囲の流体流のシミュレーションを行うと，装置の下に乱気流が形成される．その原因となるのは，[Flores-Castillo, 2022]で言及されているように，モデルの形状と流体流の状態である．流動様式サブセクションでは，より理論的に解説している．

PDCの下の乱気流の生成を避けるために，PDC装置の下に固体領域（赤）を加えた．この主な機能は，PDCの下に乱気流が生成されるのを防ぐことである[Flores-Castillo, 2022]．さらに，乱気流の生成を避けるためには以下の図に示すように，空気領域の右境界にPDCの右角を置く方法がある．

161.gif

PDC装置を灰色で，固体領域を赤で，周囲の空気領域を青で示した，修正した解析の領域Ω．

PDC装置の長さ は1，厚さ は0.0253である．PDCの厚さは，その10% が凝集器，残り が断熱材である．また，PDCは水平から30°傾いている．

168.gif

PDC装置の寸法．凝集器が白で断熱材が灰色である．

まずPCD領域を，軸に対して対角の矩形として指定する．それから幾何学変換によって，矩形は{PDCcenterX,PDCcenterZ}を中心として水平から30°回転される．

領域の長さ，厚さ，半分の厚さを指定する：
凝集器と断熱材の厚さを定義する：
PDCの長さと厚さ，その半分の長さと角度を定義する：
凝集器と断熱材の座標を定義する：
PDCの中心座標を指定する：

座標を幾何学的に変換するために，RotationTransformとTranslationTransformを使う．

変換関数を作成する：
回転し移動する：
PDCの回転と移動を可視化する：

灰色の矩形は断熱材であり，赤い点が付いた方は移動し回転したものである．同じ回転と移動を凝集器にも適用する．これは断熱材の上部にある黒く細い線として見える．凝集器は断熱材と比べると非常に小さい．

メッシュ生成

境界メッシュはToBoundaryMeshを使って生成する．PDCの座標に加え，PDC下部の乱気流の発生を防ぐために固体領域を定義する．新しい固体領域は三角形であり，PDC装置の下に置く．三角形の座標のうちの2つは断熱材の底の2つの座標であり，3つ目の座標は領域の右境界の特定の高さ にある．

最後の座標の 位置を定義する：
固体領域の最後の座標を定義する：
シミュレーション領域の境界メッシュを設定する：
点要素を持つ境界メッシュを可視化する．固体領域の中心点は赤，3番目の点は青で描画される：

PDE上部の領域は非常に重要である．その領域で，熱移動係数 の計算（方程式7）についての微分温度 が得られる．正確な結果を得るために，その特定の場所ではより細かいメッシュが使われる．その領域を微調整するために使われるメッシュ調整関数は，厳密にPDC表面の上部の半円で与えられる．

調整関数および要素によって使われる領域を定義する：

上記は，調整領域内では0.00005単位より小さい面積のメッシュ要素を，調整領域外では0.0001単位より小さい面積のメッシュ要素を生成するメッシュ調整関数を定義する．これで完全な要素メッシュが生成できる．

境界メッシュから完全メッシュを生成する：
結果のメッシュを可視化する：

材料パラメータ

流体媒体は空気であり，次にその材料特性を設定する．

空気の密度 と粘度 の値を定義する：

速度関数

このセクションでは，近くの測候所から得られた流入風速に基づいて，PDC周囲の速度場を計算する．風速は ほどである．それらを直接使うとNDSolveの実行時間が不必要に長くなる．高速度であることで，PDE，NDSolve，その他のFEMソルバの運動が増加するため，PDEを時間積分するのにより小さい時間刻みを取らなければならないからである．

これを避けるために次のアプローチを使う．ナビエストークス方程式で約 あたりの遅い流入速度および時間 までについて解く．次に， で計算した速度場以上の速度場を外挿する．これはPDC表面上の流動様式が層流であるとすると妥当である．このアプローチは流入速度が増加するにつれて速度場が線形に増加すると仮定している．このアプローチは[Flores-Castillo, 2022]で立証されており，流体モデルを追加の後処理ステップとして解いた後に説明し実行する．

流入境界条件と初期条件が一致するために，時間経過とともに流入を増加させるヘルパー関数 を生成する．平滑化したステップ関数は次のように定義できる：

ここで関数 の最小値と最大値は ，で表される．ステップの位置は で，平滑化された傾斜は で制御される．

時間経過とともに流入境界条件を増加させる関数を定義する：

ランプ関数に選ばれるパラメータは[Flores-Castillo, 2022]から取る．

ランプ関数のパラメータを定義する：
使用されるランプ関数を定義し可視化する：

次のステップでは，z軸に沿った流入風速プロファイルを生成する．速度プロファイルは方程式(8)で与えられた従来の対数変動に従う：

U(z)=U10 Log(z/zc)/Log(10/zc)

PDCは厳密に地面に接してはいないので，速度プロファイルは 方向に だけシフトする．

速度プロファイル関数を定義する：

における速度は，10 における既存の気象的風速より低い値になる．速度はより小さい流入速度について解かれる．

における速度は であると定義する：

値0.35は，シミュレーションの収束を高める効果があるため選ばれた．このモデルを使って他の値もテストされたが，値0.35が最善であった[Flores-Castillo, 2022]．

速度プロファイルを可視化する：

流動様式

レイノルズ数の値によると，流れが層流または乱気流の場合，次のように計算できる：

ここで は流体の密度， は粘度， は平均速度， は代表長さである．

流れのタイプが判別できるように，常に流体流のシミュレーションを行う前にレイノルズ数を計算することが望ましい．開放された平面構造の場合，乱気流は，つまり の風速 で発生する．これはPDC装置周辺の空気の流れのほとんどが層流ということである[Beysens, 2018]．それにもかかわらず，これは理論的概念であり，PDC形状の小さい不均衡性が流れの乱気流をよりずっと小さい値にする．こういう訳で，モデルにはシミュレーションで乱気流が発生するのを避けるためだけにPDCの下に固形領域がある．

このモデルでは，平均速度は における速度プロファイル，の 座標によって与えられ，代表長さはPDCの長さ とされる．

における速度プロファイルを使って流れのレイノルズ数を評価する：

この値により，理論的にこの流れは層流であることが分かる．

境界条件

流体力学モデルには4種類の境界条件がある：

  • 流入
  • 流出
  • 壁面境界条件
  • 上方境界条件

入り口では流速はに設定される．は速度プロファイルである．

入口境界における流速を設定する：

領域の地表レベルは壁面およびPDC装置の外部境界として表されるので，滑りなし流れ状態をモデル化するために，流速は境界上ではゼロに設定される．

PDC装置の壁面で指定される滑りなし条件は，表面の速度が0である必要があるという条件を満たすために使われる．境界における0速度は，凝縮モデルのセクションで説明してあるように，熱移動係数 の計算には必須である．

領域の下部と壁面境界で滑りなし境界条件を設定する：

出口では圧力出口境界条件を使って，流出する流れをモデル化する．ここで出口圧力は絶対圧力，に設定される．これは相対圧力 に等しい．

流れ出口における圧力出口境界条件を設定する：

領域の上部では，開境界条件が適用される．開境界条件は幾何形状の流体部分が終るが形状で表されていない残りの流体領域に接している境界のことである．この境界条件はノイマン0値としてモデル化される．これらは自然なデフォルト境界条件なので，考慮しない．

境界条件を組み合せる：

PDEモデルを解く

まず，流体力学PDEを流速場 について解く．

2Dの流体力学モデルを構築する：
初期条件を定義する：
モデルパラメータを含む流体力学PDEを定義する：
シミュレーション時間を設定する：

ナビエストークス方程式は高指数微分代数方程式なので，方程式を効率的に時間積分するオプションが指定される．時間積分器の最大微分次数は，拒絶されるステップ数を最小限にするために還元される．時間依存境界条件は，時間積分するのにより効率的な一貫した方程式系を作成するために微分される．最後に，速度が圧力より高い次数で補間されると，安定解を見付けることができる．ここで速度 と は二次で補間され，圧力 は一次で補間されるよう設定される．

有限要素についてNDSolveで使えるオプションおよびそれらを指定する方法については，チュートリアル有限要素のためのNDSolveオプションに記載されている．

また線形ソルバとして，NDSolveがデフォルトのアルゴリズムではなくマルチフロンタル法を使うよう指定した．この場合，デフォルトの線形ソルバは与えられたほぼ特異行列を解くことができない．デフォルトを使うと，次のメッセージが現れる：

流体力学PDEモデルを解く：

流体場を調べるために，以下の可視化で速度場 の流線プロットおよび速度の大きさの等高線プロットを示す．

後で行う可視化のために固体領域を構築する：
における速度場の流線と大きさを可視化する：

速度場の外挿は次のサブセクションで行う．

外挿法

ナビエストークス方程式は約の低い流入速度について解かれ，時間まで解いた．しかし，熱移動シミュレーションを行った全体の時間，つまり時間またはをカバーすることができ，実際の風速値を持つ解が欲しい．これを達成するために，前に得た風速場を補間して，すべての時間で実際の風速にマッチさせる．いつ速度値を修正する必要があるかを知るために，すでにここで得た風速関数を使う．

補間がどのように動作するかを理解するために，仮想の の場合を見てみる．例えば地表からの風速 がであり，これをでの気象的速度場の値と比較する．この場合のように，もし風速が より大きければ，におけるナビエストークス解からの解の場は，因数で乗算され，新しい速度場はの場になる．速度が より小さい場合は，の速度場は指定された におけるナビエストークス開から直接取られる．この時間刻みにおけるナビエストークス方程式からの速度場はにおける風速の大きさをすでに再構築しているため，ここでは での解を取ることにする．

解の場は まで解かれるので，速度が より小さいときは，より大きい時間について，ナビエストークス方程式からの解がまだ利用できる．この時間の不整合性を修正するために，2つの関数を構築する．

最初の関数は入力として と風速データを取り，出力として流体シミュレーションの範囲の間にある時間を与える．この関数は，例えば等の長い時間のときにどの時間速度解を使えばよいかを知ることを目的としている．この関数によるとにおける速度解が使える．

取得する時間は以下の方程式を解いて得る：

ここで は流体シミュレーションで使ったのと同じランプ関数であり， は定数値である．

ランプ関数は について解く関数であり，入力として与えられる値 はまたは割合である．値は の値，そして結果的に時間の戻り値を決定する閾値でもある．

Solveを使って方程式を記号的に解く：

Solveは常にネストされた規則のリストを返す．解を使うためには，まず規則を置換してリストから値を抽出する必要がある．

この関数は割合とを比較する．速度の割合が0.98以上ならば，とするとランプ関数は について解く．条件が満たされない場合，とするとランプ関数は について解く．

条件が真なので，戻り時間値は となる．それ以外の場合，速度の割合に依存する低い時間値となる．

使用する時間値を割り当てる関数を作成する：

この関数を定義したら，流体シミュレーションと熱移動シミュレーションの時間を関連付ける補間関数を作成する． 値は0からの範囲， 値はからの範囲である．この関数を使うと，熱移動解の時間領域にかかわらず，より長い時間についての解を抽出することができる．

時間補間関数を定義する：

すでにナビエストークス方程式の解にアクセスする方法があるので，補間の実装に進むことができる．

実装ではIf条件文を使って，風速が より大きいときと より小さいときを比較する．時間 における風速が ()における速度より遅ければ，時間 における速度場は修正なしでナビエストークス方程式からのものになる．そうでない場合は，時間における速度場は，時間 における風速と における速度の割合で乗算される．

ここで の範囲はからの間であり，time はとの間の時間値を返す関数 timeFunction である．

方程式17からの速度 は以下のようになる：

において補間された速度場を可視化する：

最後に，速度関数を定義するために，熱移動モデルはPDC領域を含む異なるメッシュを使うということを考慮しなければならない．したがって，PDC装置は固体であるためそれに割り当てられた速度場がない．このため，後に定義する要素マーカーを次の関数で使って，PDC内部で速度ゼロを割り当てる．

熱移動モデルで使われる速度を作成する：

熱移動モデルの設定

概要

このセクションでは，方程式，解析の領域，材料パラメータ，境界条件を含む熱移動モデルを説明して設定する．

熱方程式は流体領域でもPDC装置内でもアクティブである．過渡熱移動モデルの場合，温度分布 は過渡熱方程式で記述される．過渡熱方程式は以下で与えられる：

モデル変数を設定する：

この次は一般的な熱方程式であり，シミュレーションで使われるものではない．パラメータは熱移動パラメータの指定サブセクションで指定される．

一般的な熱方程式を計算する：

熱移動PDEモデルの他に，境界条件を指定することも必要である．

  • 対流項 は装置上の空気流れ速度をモデル化し，PDCから熱エネルギーを除去する．
  • 放射 は放射を介した周囲とのPDCの熱エネルギ―の流入と流出をモデル化する．
  • 流入熱流束と考えられる凝縮過程 も装置の温度に影響を与える．
  • 流入熱流束と考えられる太陽放射照度 も装置の温度に影響を与える．

要するに，熱移動モデルはさまざまな境界条件を考慮に入れて，流体の温度のシミュレーションを行う一方で凝縮器と断熱材の温度分布のシミュレーションも行う．

注意すべきことは，この設定はPDE装置自体のサポート構造は考えていないということである．凝縮表面は断熱材に固定されている[Beysens, 2018]ため，マウントのシミュレーションは不要なのである．

これ以降は，シミュレーションすべてで気象データを使う．

メッシュ生成

このサブセクションでは，メッシュ要素のサイズを増大することによってメッシュ要素の数を削減する．このメッシュには流体シミュレーションでは考慮されないPDC領域が含まれた．流体シミュレーションと同様に，赤い固体の三角形領域も削除された．この領域はシミュレーションの関心領域ではないので，ここを削除しても熱方程式の解は変わらない[Flores-Castillo, 2022]．

境界要素メッシュを生成する：
境界要素メッシュを，その点要素と中心点を使って可視化する：

前述の通り，この領域はそれぞれ異なる特性を持つ材料からできている．そのため，異なる部分領域にマーカーを割り当てる．これらのマーカーは材料パラメータを設定するために使われる．これにマーカーを使う方が，各部分領域に式を指定するよりも簡単である．

マーカーが異なる部分領域に割り当てられるように，これらの部分領域にある座標を指定する必要がある．

マーカーの座標を設定する：

要素マーカーは"Air"->10，"Condenser"->20，"Insulator"->30である．

要素マーカーを割り当てる：

異なる部分領域を指定するとき，それぞれの領域を微調整することもできる．

各部分領域に対するセルサイズを定義する：
完全メッシュを生成する：
新しいメッシュを可視化する：

上記のメッシュから分かるように，メッシュの調整はPDCの表面と空気の間の接触面だけに集中している．

空気と凝縮器の間の境界を調べる：

青い領域は空気，灰色は断熱材，白は凝縮器を表している．凝縮器の表面を調整する必要はない．関心があるのは凝縮器表面と空気およびその上の領域の間の接触面である．

熱移動パラメータの指定

各部分領域の材料特性を指定する必要がある．流体媒体は空気である．凝縮器は低密度ポリエチレン(LDPE)，断熱材は発泡スチロールである．

各領域について比熱容量 を設定する：
熱伝導率 を設定する：
各材料について密度 の値を指定する：
熱移動モデルのパラメータを設定する：

それから外挿法サブセクションで説明た速度を指定する．

"HeatConvectionVelocity"エントリを指定する速度を設定する：

境界条件

熱移動PDEモデルの他，以下の境界条件を指定する必要がある．

  • 対流項 は装置上の空気流れ速度をモデル化し，PDCから熱エネルギーを除去する．
  • 放射 は放射を介した周囲とのPDCの熱エネルギ―の流入と流出をモデル化する．
  • 流入熱流束と考えられる凝縮過程 も装置の温度に影響を与える．
  • 流入熱流束と考えられる太陽放射照度 も装置の温度に影響を与える．

次に熱移動モデルの境界条件項を数量化する．周囲は ，凝縮器は ，空は という略称を使う．これらの境界条件は内側境界に適用されるものであり，シミュレーション領域の境界に適用されるものではないという点に注意が必要である．

対流項

PDC装置の壁で滑りなし流体条件を指定したため，表面の速度は0である．したがって，対流境界条件は，熱移動モデルで直接使うことができない．それでも，凝縮モデルセクションで説明するように水の量を計算するのに必要な熱移動係数 は計算する．このプロセスは対流が考慮に入れられることを確実にする．

凝縮項

シリアルアプローチでは，PDCの温度は凝縮項 を使わないので，この境界条件は使われない．

放射項

は放射冷却項であり，天空から吸収された放射と凝縮器によって発せられた放射の間の差分として表される．これは以下の方程式で表せる[Flores-Castillo, 2022]：

ここで は天空の放射率， は凝縮器の放射率，はシュテファン・ボルツマン定数， は範囲が0から1の形態係数である．

方程式(9)に従って天空温度 が周囲温度 に関連していると考えることにより， はさらに簡約することができる：

Tsky=Taϵsky1/4

方程式(10)を(11)に代入すると以下が得られる：

[Flores-Castillo, 2022]で与えられている方程式(12)で表される境界条件は，通常使われるHeatRadiationValueにマップしないという意味で従来型ではない．それでも与えられた境界条件は，実装することができる．ここで示される放射冷却モデルは，通常入射を表すときに使われる凝縮器表面の吸収率 α を考慮していない．詳細は熱放射境界条件または[Beysens, 2018, equation 7.20]でご覧いただきたい．ここではHeatRadiationValueを使ってモデル化することができるより一般的な形式の方程式が示されている．

このサブセクションでは，放射項 をさらに調整する．上の方程式を使う代りに，天空温度 に対する雲量比の影響も考え，を にする．

この境界では，雲量 ，放射率 ，シュテファン・ボルツマン定数 ，PDCの傾斜角度に依存する形態係数 という変数を使う．

形態係数

形態係数 の値を得るために，スペクトル放射輝度 を記述し，さまざまな波長 において表面要素によって反射される力を提供するプランクの法則を使う．プランクの法則をは物体の単位表面積あたり，単位立体角 あたり，単位波長あたりで測定される[Beysens, 2018]．プランクの放射法則は以下である：

さらに計算するためにPlanckRadiationLawを使う．PlanckRadiationLawは方程式(13)と同じ公式を使う．

スペクトル放射輝度 が，与えられた温度のすべての波長および立体角全体で積分されると，完全な黒体の放射輝度が得られる：

ここで半球上で を積分し，からで を積分する．黒体はランベルトの余弦則に従うので，余弦が現れる．この方程式では は天頂角()を， は方位角()を表す．

を含む積分の部分はPlanckRadiationLawで計算され，残りの部分は記号積分を使って手動で計算される．

温度 について積分を計算する：

PDCが傾くと，立体角は全天空上では積分されず，計算は2つの積分に分割される．最初の積分 は0から までと から までの範囲である．2つ目の積分で変更されるのは の範囲であり，傾斜 θ から である：

上の方程式を計算する：

指定された角度に傾いたPDCの形態係数は と の比で与えられる．

形態係数 を計算する：
シュテファン・ボルツマン定数 σ
シュテファン・ボルツマン定数 を定義する：
凝縮器の放射率

前述の通り，放射冷却は水を集めることを可能にする必須の熱工程であり，装置の放射率 はPDCからどれだけのエネルギーが反射されるかについて重要は役割を果たす．表面の放射率は温度，波長，角度，その他の変数等いくつかの要因に依存する．ここではモデルを簡単にするために一定の放射率を使う．

冷却量の値 は凝縮器の放射率 の値に敏感であることを知っておくことは大切である．したがって，対照実験のデータにフィットするようにシミュレーションの中で放射率の値を調整することができる[Beysens, 2018]．

このモデルでは凝縮器の状態によって2つの異なる放射率を使う．凝縮器に水がある場合は水の放射率 を使い，それ以外の場合は凝縮器の放射率 を使う．ここでは両方の放射率について同じ値を使っている．これらの値は常に0と1の間でなければならない．

LDPEと水の凝縮器表面の放射率の値を設定する：

の定義の中でIf条件文を使い，表面温度が露点をいつ下回るかを知り， か かを選ぶ．露が表面に形成されるためには，表面温度は露点温度 より低くなければならない．

を設定する：
境界位置

太陽放射照度の熱流束 と熱放射 はともに凝縮器表面で起こる．

凝縮器表面は線 を使って指定する．線 は表面の 座標を与える．凝縮器表面に重なるこの線を生成するためには，線の方程式 を使う関数を定義する．ここで は線の傾き， は 軸との交点である．

線の方程式を使って関数を定義する：

の入力を以下で定義する．

線 と入力角の 交点を定義する：
表面境界を定義する：
境界条件を適用する境界を可視化する：

は，HeatRadiationValueが使うものに直接マップしないので，熱放射境界条件はHeatRadiationValueで実装することができない．このような訳でNeumannValueを設定する．

NeumannValueを使って熱境界条件を設定する：

太陽放射照度

凝集器に熱が加わる別の要因として太陽放射照度 がある．太陽放射照度は，日中にしたがって凝集器の温度を上げたり下げたりする．このモデルでは簡単な太陽反射モデル[Wilkerson et al., 2013]を使う．この太陽放射モデルでは，地球上のあらゆる場所の日の出日の入り時刻を考慮して，太陽からどれだけの熱流束が来るかを予測することができる．これは後でノイマン境界条件として実装する．

まず回転行列として表される，地球の回転軸を定義する．ここで回転角は地球の傾斜角()である．

傾斜角をラジアンで定義して，軸回転行列を構築する：

以下の2つの方程式は地球の回転を与える．最初の方程式はで与えられる回転角を持つ回転行列であり，この角は季節や時間によって異なる．

回転角と地球の回転行列を定義する：

合計時間は秒で与えられ，正午の0から始まる[Flores-Castillo, 2022]．

季節は1月1日から最終日までの日数をで割ったものである．

指定された日付に従って季節を定義する：

1年を通して太陽に対する地球の位置は以下の関数で定義する．

地球の位置ベクトルを定義する：

次に地球の位置をベクトルとして定義する．

位置ベクトルを定義する：

2つの回転行列と位置ベクトルから，表面に垂直なベクトルを得ることができる．

法線ベクトルを定義する：

最後に，入射角に最大太陽放射照度1360 を掛けて，入射放射を計算する．

入射放射関数を定義する：

太陽放射照度は，利用可能な実験データにフィットするように調整された．この調整は最初の時間に行われたので，シミュレーションの日の出・日の入り時刻は実際のものと合致する[Flores-Castillo, 2022]．

最初の時間を調整する：
位置から緯度を得る：
雲量と凝集器の傾斜角を考慮して，指定された日時におけるミルレフトの太陽放射照度を可視化する：

太陽放射照度のプロットは，日の入りと日の出がそれぞれおおよそとに起きていることを示している．太陽放射照度は夜はゼロになる．

期間内に蓄積された太陽放射照度を可視化するために，太陽放射照度関数から点のリストを作成する必要がある．その後で関数Accumulateを使って合計の太陽放射照度を得る．

熱流速関数 を設定する：
点のリストを作成する：
Accumulateを使う：
この期間内に蓄積された太陽放射照度を可視化する：

熱流束境界条件は太陽からPDC表面に流入する熱エネルギー量を表す．

PDC表面における境界条件を設定する：

ここで"HeatConvectionVelocity"を削除したのは，この境界条件が固体のPDC表面のためのものだからである．熱対流速度が含まれていれば，HeatFluxValueはここでは必要ない境界単位法線 に関連付けられた対流速度項を含む流体の境界をモデル化することになる．

HeatFluxValueに関連付けられた対流速度を持つHeatFluxValueの出力を表示する：

入り口温度

これは，シミュレーション領域境界に適用される唯一の境界条件である．

流れの入り口 において，温度は周囲温度 に設定する．

入り口における温度表面境界条件を設定する：

解と可視化

このサブセクションでは熱方程式を解く．まずすべての領域上で周囲温度 を表す初期条件を定義してから，PDEを設定する．

初期条件を定義する：
モデルパラメータで熱移動PDEを定義する：
シミュレーション時間を設定する：

また，このモデルでは方程式を時間積分するためのオプションが指定される．MaxStepSizeは900，StartingStepSizeは1が選ばれた．使われる時間積分法は，"TimeIntegration"オプションで設定され，最大微分階数が1の"IDA"メソッドである．最後に時間依存境界条件を微分して，時間積分を行うのにより効率的な一貫した方程式系を作成する．NDSolveが有限要素について持つオプションおよびその指定方法については，有限要素のためのNDSolveオプションチュートリアルに記載されている．

熱移動PDEモデルを解く：

この警告メッセージは，使用されたモデルとメッシュの性質のため，生成されている．これはモデルの結果には影響しない．

における温度分布を可視化する：
における凝集器表面付近の温度分布を可視化する：

上のプロットから，凝縮器表面の温度と凝縮器の上の空気温度が分かる．どちらの温度も熱移動係数 および凝縮した水の量 を計算するために必須である．

凝縮モデル

概要

凝縮モデルは過飽和方程式(14)に基づいている：

この方程式には凝縮と蒸発の2つの過程がある．の値が正ならば凝縮が起こり，となる．の値が負ならば蒸発が起こる．このモデルでは蒸発を考えないので，質量変化率はゼロ と設定する．

ここでゼロとは水の蒸発は考慮しないということである[Flores-Castillo, 2022]．したがって，凝縮した水はすべて集められるので方程式は以下のように簡約できる：

方程式(15)は以下の積分として見ることができる：

この積分では凝縮した水を計算することができる．熱移動係数 を除いて，方程式(16)の値は後で示す公式から直接得ることができる．この段階では熱移動係数 の計算方法をもう少し説明する必要がある．

熱移動係数 の計算は，PDC上の流れが層流であるということと表面の流速が0であるということの2つの仮定に基づいている．

流体における熱方程式の導出で説明しているように，伝導成分 と対流成分 がある．対流成分は流速 に比例する．流体とPDC表面の間の接触面における熱移動係数 が求めたいが，境界における流体流速を0と設定したので，流速の対流部分は境界においてゼロでなければならない．PDC表面では，壁にフーリエ熱伝導法則 を適用したままにする．

境界ではフーリエ熱伝導法則をニュートンの冷却の法則と同じであると考えることができる：

ここで は熱移動係数，は周囲温度，は凝縮器温度である．

と見なすことによって，熱移動係数 が与えられる：

流体の流れが層流ではない，または流体の流速が境界においてゼロではないならば，熱移動係数 の計算には異なるモデルを使う必要がある．

熱移動係数 と凝縮器の温度 を知ると，凝縮速度 を計算することができる．

このモデルは凝集した水がすべて収集され，水の蒸発はないと仮定していることをここでもう一度強調しておく．

この連続アプローチは凝縮した水の量を推定することはできるが，水の凝縮がある場合PDCの表面温度 は低くなるため，高く見積り過ぎる可能性がある[Flores-Castillo, 2022]．

潜熱なしの水生成

凝縮された水は，まず方程式(17)を使って熱移動係数を計算し，次に方程式(18)を積分することで計算できる．

熱移動係数

熱移動係数を計算するためには，まずPDCに垂直な温度勾配 を計算する必要がある．

温度勾配はPDC表面と表面に平行な線の間の温度および距離の差を得ることで近似される．

勾配を計算するためには，PDC表面を200のセクションに分け，平面において解から直接温度を取る．2つ目の温度はPDC表面に平行な線上の距離から取る surface is divided into 200 sections, and the temperature is taken at the plate directly from the solution. The second temperature is taken from a distance over a line parallel to the PDC surface [Flores-Castillo, 2022].

短い距離を定義する：
セクションの数を定義する：

平行線を作成し，その戦場で評価するために，直線方程式および位置xに依存する 軸との新しい交点を作成する．

線の新しい方程式を定義する：
PDC表面に平行な線の 交点を定義する：

次に，表面に垂直な平行線で温度評価を行う関数と，温度差を得る2つ目の関数を定義する．

両方の関数を定義する：

最終ステップは，時間とPDC位置の関数として の補間関数を生成することである．

補間関数として を計算する関数を定義する：

と についても同じことを行う．

補間関数として を計算する関数を定義する：
補間関数として を計算する関数を定義する：

この関数はシミュレーションの指定された時間または一定期間で使える．その後補間関数を使って熱移動係数とその他の変数の挙動を可視化することができる．

まず，における熱移動係数と温度勾配の挙動を示す．

における 関数を設定する：
熱移動係数の補間関数は空間だけに依存することを確認する：
平均の熱移動係数を計算する：
における 関数を設定する：
での熱移動係数とPDC表面 の温度差を可視化する：

このプロットから，関数が粗いことが分かる．これはメッシュを細かくしたりデータをフィルタすることで向上させることができる．後者はGaussianFilterを使って適用する．

の終始の挙動が得たい場合は，それぞれの時間を表す異なる補間関数を作成する必要がある．この場合， の計算は毎に行われ，これらの補間関数から時間と空間に依存する一つの補間関数を作成する．

時間刻みを定義する：

補間関数の時間領域はからまでである．ゼロによる割り算を避けるために時間は省略する．

時間と位置に依存する の補間関数を作成する：
温度勾配 と表面温度 についても同じことを行う：
における平均の熱移動係数を計算する：
ガウスフィルタを使ってにおける熱移動係数 を可視化する：

このプロットから，空気が最初に触れる左の表面の方が高い熱移動係数を持つことが分かる．このプロットは，局所化された熱移動係数に基づいてPDC表面のある部分がより多い水を凝縮するかを示している[Flores-Castillo, 2022]．

ガウスフィルタを使って における温度勾配 を可視化する：

計算した熱移動係数が有効であるかどうかを知るために，その熱移動係数と理論的近似を比較することができる．それが以下の方程式である：

ここで は補正係数， は実験値， は典型的な凝縮器の長さスケールである．

速度 板上および板付近の風速である[Beysens, 2018]．したがって，地上1mの風速を使う．

=1，=4，の値で熱移動係数を定義する：
熱移動係数近似を可視化する：

熱移動係数の範囲は0から5であることが分かる．

水の凝縮

水の凝縮は方程式(19)で推定することができる．方程式を解くのに必要な圧力は両方とも以下の式で計算する：

ここで は蒸気圧を表し，（周囲温度 または凝縮器温度Tc [°C]）次第で または になる．

蒸気圧を定義する：

最後の未知数である蒸発の潜熱 を計算するために，温度範囲248 から313 の間に適用される以下の方程式を使う：

温度範囲が273 から290 の間であれば，方程式(26)は簡約できる：

この計算には方程式(20)を使う．入力温度はで与え，出力温度はで与える．

蒸発の潜熱を定義する：

これらの関数を使って，方程式(21)に従って の計算に進むことができる：

さらに，ここで結露に必要な温度条件を加える．それは「表面温度 が露点 以上ならば，質量速度はゼロになる」というものである．

湿度定数を定義する：
質量速度方程式を定義する：

表面積 は1であり，蒸発の潜熱 に1000を掛けてSI単位を得る．

さまざまな変数がどのように振る舞うかをよりよく理解するために，それらを異なる時間で可視化することができる．

における温度 ，露点 ，周囲温度 を可視化する：

表面温度のプロットを見ると，では が露点 より低いため，凝縮器表面の において水蓄積が起こることが分かる．

における水蓄積率を可視化する：

このプロットはまさに で水が蓄積され始めることを示している．

最も一般的な降雨量測定は，指定された期間における総雨量であり，容器内の水が到達した最高点を通常ミリメートルで表す．どれだけの体積の水が凝集したかを知るためには，集水域を知る必要がある．集水域が1ならば，1は1に等しい．凝集率はまたはで測定される．

方程式(22)を積分するために，

を計算する．以下で説明するように近似を適用する．

まず，NIntegrateを使って時間 においてPDCの長さに沿って方程式を数値的に積分する．次にこの結果に使用した時間刻みを掛ける．各時間間隔について値のリストを生成し，最終的にAccumulateを使って凝集した水の総量を求める．

の積分の時間刻みを定義する：

最初の時間刻みでは，NIntegrateは収束しない．この時間の の値が表面全体でゼロの値だからである．NIntegrateはエラーメッセージを出力し，これらの時間刻みではゼロの値が返る．

次の計算では，NIntegrateが生成するエラーメッセージは表示されない．

値のリストを生成する：
からまでの時間のリストを作成する：
の形式で値を入れる：
補間関数を作成する：
凝縮した水を可視化する：

このプロットから，凝縮の潜熱が存在しない場合の計算された水の凝縮は約0.226，つまり226に等しい．

シミュレーションでは凝縮項 は使われなかったが，方程式に応じて計算することができる：

ここで は凝縮の潜熱， は凝縮水の質量流量である．凝縮の潜熱は水が気体から液体に変化するときに放出される熱である．

はすでに計算してあるので，残りは凝縮の潜熱 だけである．これは次の方程式を使って計算できる：

ここで凝縮の潜熱 は蒸発の潜熱 の負の値である．

凝縮の潜熱 を定義する：
を関数として定義する：
における凝縮の潜熱によって放出された力を可視化する：

このプロットは， 凝縮器表面の水の凝縮によって，において最大値おおよそで熱が凝縮器表面で効率的に放出されることを示している．時間が異なると異なる量の熱が放出される．

検証

[Lekouch et al., 2012]によると，モロッコのミルレフトで観察された水の凝縮は0.22，つまり220であった．対照的に，凝縮の潜熱がない計算された水の凝縮はほぼ0.226，つまり226であった．水量を比較する他に，PDCの表面温度と[Lekouch et al., 2012]の実験結果を比較することができる．

実験データから補間関数を構築する：

PDC表面では，表面温度は について異なる．したがって，表面の左，右，中央という特定の場所における温度を表す3つの異なる表面温度を示す．また，平均の表面温度はプロットに含まれている．

異なる温度を可視化する：

このプロットから，表面温度 が急激に上昇および下降するときのプロットの2つの極値において太陽放射照度の影響が観察できる．これらの高い表面温度は太陽放射照度によるものなので，夜にだけ水がPDC表面に凝縮できるということが分かる．

からまでの範囲で異なる温度を可視化する：

このグラフが示すもう一つの特徴は，表面の左側の温度が右側よりも常に高いということである．これは熱移動係数 の値が表面の左側の方が高いためである．また，雲量の影響も 付近で観察できる．雲量の影響が大きくなるため表面温度 も上がるのである．

の後の領域の表面温度の値は観察された実験値により近くなる[Flores-Castillo, 2022]．

付録

WeatherDataは世界中の気象データへのアクセスを可能にする組込み関数である．ある期間が与えられると，この関数は気象パラメータの時系列を返す．測候所がない場所では，任意のモデルの周囲条件の妥当な推定として，近隣の測候所のデータを使うことができる．

異なる測候所からのデータへのアクセスを可能にする関数を定義する．WeatherDataは時系列として戻るので，これらのデータ点を連続関数，つまり補間関数に変換することが必要である．欠損データの場合，MissingDataMethodを指定する必要がある．

この場合，MissingDataMethodは欠損したデータ値を埋めるために，0次の補間を使う．その一方で，補間次数2でInterpolationを使うことによって補間関数が生成される．

補間関数と最後のデータの秒数を返す関数を定義する：
ある場所の雲量率の気象データを見る：
関数を適用して補間関数を生成する：

測候所からの気象データ

ここに挙げるのは関数 WeatherParameter がどのように作用するかを示す例である．使用される場所はモロッコのミルレフトに一番近い測候所である．この場所は，標準の測候所の識別子であるICAOコードとして与える．

場所を定義する：
緯度と経度を抽出する：
測地位置を定義する：
高度をとで得る：
日付をDateオブジェクトとして定義する：
WeatherDataの期間を定義する：
両日の間の経過時間をで得る：
GeoGraphicsで場所を可視化する：
雲量，相対湿度，風速，周囲温度，露点に対する補間関数を作成する：

WeatherDataには，定義した期間のデータがない場合がある．例えば，ここではもとの23.75時間の代りに21時間分しか得られない．

パラメータが同じ期間をカバーしない場合もある．これを修正するためには，データの最小の長さを定義する．変数のプロットや他の計算を行うときにこれが使われる．

すべての気象パラメータの最小の長さを定義する：
周囲温度 ，上空温度 ，露点 ，雲量の影響のある上空温度 をプロットする：
雲の影響のある上空放射率 と雲の影響のない上空放射率 の差分を可視化する：
相対湿度と雲量比を可視化する：
との風速をプロットする：

用語集

参考文献

1.  Beysens, D. (2018). Dew Water. Aalborg: River Publishers.

2.  Flores-Castillo, Pedro. (2022). Computational Analysis of Water Condensation for Passive Dew Condensers to Optimize Water Condensation. Master's thesis, Harvard University Division of Continuing Education.

3.  Gires, A. (2018) How Do We Measure Rainfall?. Front. Young Minds. 6:38. doi: 10.3389/frym.2018.00038

4.  Lekouch, I., Lekouch, K., Muselli, M., Mongruel, A., Kabbachi,B., & Beysens, D. (2012). Rooftop dew, fog and rain collection in southwest Morocco and predictive dew modeling using neural networks. Journal of Hydrology (Amsterdam), 448-449, 60-72. doi: 10.1016j.jhydrol.2012.04.004

5.  Wilkerson, S., Wattenberg, F. A., & Park, R. (2013). Solar Energy Incident on Earth's Surface. Wolfram Demonstrations Project. Retrieved from https://demonstrations.wolfram.com/SolarEnergyIncidentOnEarthsSurface/