\(\require{physics}\)

Outflow境界の設定方法

例:貯水池の放流

貯水池の放流には表層放流や選択取水をはじめ,比較的複雑なオペレーションが利用されることが多い.これらを実現する方法を説明する.利用するオブジェクトはPointFluxBoundary,もしくはColumnFluxBoundaryとなる.
PointFluxBoundaryは,指定した計算セル(Box)において流入出(流出:fluxはマイナス,流入:fluxはプラス)させるオブジェクトであり,与えた時系列ファイルに従って動作する.個々のセルを任意の時系列で操作可能であるため,(面倒かどうかは別として)多くの状況を再現できる.一方のColumnFluxBoundaryはコラム単位で流出させる(Columnに存在するセルの体積の割合で流量を分配)ため,鉛直方向に分布を指定して流入出を再現できない(流出:fluxはマイナス,流入:fluxはプラス).しかしながら,下限(lower_limit)を指定することで,表層放流を近似的に表現できる.これらのオブジェクトは(セル,もしくはコラムが)干出する(dryになる)場合に流入出が止まることに注意する必要がある.データの与え方は河川流入などに利用するColumnFluxBoundaryについてと同じであるが,ColumnFluxBoundaryの場合はupper_limitとlower_limit(もしくはいずれか)をfluxと同じように時系列で指定することで流出させる高さ・範囲を変更することができる.upper_limitを省略した場合はupper_limitは水面,lower_limitを省略した場合はlower_limitは水底となる.ただし,このColumnFluxBoundaryの方法はセルの体積に比例させて流量をセルから抜き取っているため,近似的な流出方法であることに注意(境界に流速を指定しているわけではない).

時系列データの例(タブ区切り,一行目の項目名は必要,はじめの3600秒かけて線形的に増加させる場合)

time	flux	lower_limit
0	0	23
3600	-5	23
7200	-5	23
10800	-5	23
14400	-5	23
18000	-5	23
21600	-5	23
25200	-5	23
28800	-5	23
32400	-5	23

ColumnFluxBoundaryオブジェクトの生成方法は,河川流入などに利用するColumnFluxBoundaryについてを参照してもらえればよい.PointFluxBoundaryの場合は,コンテナ(global_i, global_j)とコンテナ内のcolumn(local_i, local_j)の指定に加えて,セルを指定するためのlevelであるk(k=0がdomainの一番下であり,必ずしも水底ではないことに注意する.

PointFluxBoundaryの場合

point_outflow = fantom.PointFluxBoundary(water_body, "time_series.txt", global_i, global_j, local_i, local_j, k)

ColumnFluxBoundaryの場合

surface_outflow = fantom.ColumnFluxBoundary(water_body, "time_series.txt", global_i, global_j, local_i, local_j)
surface_outflow:add_column(global_i, global_j) -- コンテナに含まれる全てのコラム
surface_outflow:add_column(global_i, global_j, local_i, local_j) -- コンテナ内の特定のコラム

また,時間進行のループ内において,陰的スキームが始まる直前に以下のように境界条件を更新する必要がある.

for counter = 0, total_count do
  time = dt*counter
  -- other explicit terms
  point_outflow:update(time)
  surface_outflow:update(time)
  -- implicit terms
end

サンプルコード

require("fantom")
require("writer")

-- utility function
function string_vector(tbl)
  local tmp = fantom.VectorString()
  for i,v in ipairs(tbl) do tmp:push_back(v) end
  return tmp
end

-- basic parameters
case_name              = "outflow_test"
dt                     = 20
total_count            = 5*3600/dt -- 5 hours
save_interval          = 600/dt    -- every 10 minutes
dx                     = 10.0
dy                     = 20.0
nx                     = 30
nz                     = 30
limit_depth            = 0.01
fantom.Parameter_THETA = 1.0

-- vertical grid
dz = fantom.VectorDouble(nz)
for i=0, dz:size()-1 do
  dz[i]=1.0
end

-- horizontal grid (vertical 2D)
fg  = io.open("geom.txt", "w")
fxr = io.open("xrfm.txt", "w")
fyr = io.open("yrfm.txt", "w")
for j = 0, 0 do
  for i = 0, nx-2 do
    fg:write ("0\t") 
    fxr:write("1\t")
    fyr:write("1\t")
  end
  fg:write ("0\n")
  fxr:write("1\n")
  fyr:write("1\n")
end
fg:close()
fxr:close()
fyr:close()

-- object generation
water_body   = fantom.WaterBody(dx, dy, dz, "geom.txt", "xrfm.txt", "yrfm.txt", limit_depth)
surface      = fantom.SurfaceSolver()
m_advection  = fantom.UltimateQuickest(string_vector({"u", "v"}))
vtk_writer   = writer.VTKWriter(water_body, case_name)
txt_writer   = writer.TextWriter(water_body, case_name)

-- object settings
water_body:set_uniform_surface_level(28.0)
surface_outflow = fantom.ColumnFluxBoundary(water_body, "surface_outflow.inp")
surface_outflow:add_column(0,0)
point_outflow   = fantom.PointFluxBoundary(water_body, "point_outflow.inp", nx-1, 0, 0, 0, 9)

-- variables
water_body:set_var("u",            0.0)
water_body:set_var("v",            0.0)
water_body:set_var("w",            0.0)
water_body:set_var("temperature",  0.0)
water_body:set_var("salinity",     0.0)

-- time marching
time = 0.0
for counter = 0, total_count do
  print("----- count = ", counter)
  water_body:construct_box_variables()
  if counter % save_interval == 0 then
    vtk_writer:write(time, counter)
    txt_writer:write(time, counter)
  end
  m_advection:apply(water_body, dt)
  water_body :update_explicit_terms(dt)
  ----- outflows
  surface_outflow:update(time)  
  point_outflow  :update(time)  
  -----
  surface    :apply(water_body, dt)
  water_body :hydrostatic_continuity()
  water_body :update_wet_components()
  time = time + dt
end

point_outflow.inp

time	flux
0	0
3600	-0.3
7200	-0.3
10800	-0.3
14400	-0.3
18000	-0.3
21600	-0.3

surface_outflow.inp

time	flux	lower_limit
0	0	20
3600	-1	20
7200	-1	20
10800	-1	20
14400	-1	20
18000	-1	20
21600	-1	20