مش بندی مدل آب های زیرزمینی با MODFLOW 6 و Flopy
MODFLOW 6 آخرین نسخه کد مدل سازی زیرزمینی MODFLOW توسط USGS است. این نسخه شامل برخی از ابزارهای نوآورانه و یک شکل مجدد کامل ساخت مدل است، اما تا به امروز این کد هیچ رابط کاربری گرافیکی (GUI) تجاری / باز نداشته است که باعث بهبود استفاده از کد برای طراحان مبتدی و مبتکران علوم زمین می شود.
یکی از جدیدترین ویژگی های جدید از MODFLOW 6 گزینه های مختلف اختیاری برای نسل مش مدل است. تنظیمات از شبکه منظم (مشابه MODFLOW 2005)، شبکه مشبک سه بعدی و شبکه غیر ساختاری است. Flopy که یک کتابخانه پایتون برای ساخت و شبیه سازی MODFLOW 6 است و دیگر مدل ها دارای ابزار برای تولید مش های سه گوش است. جریان کار در مدلسازی آبهای زیرزمینی با MODFLOW 6 و Flopy برای مدلهای مش مثلثی بسیار منعطف است و پتانسیل زیادی برای مدل سازی جریان آب های زیرزمینی محلی و منطقه ای وجود دارد.
این آموزش فرایند کامل را برای ایجاد یک مشی مثلثی با خدماتی از Flopy نشان می دهد و آن را به مدل MODFLOW 6 اضافه می کند. مدل شبیه سازی شده و نتایج به صورت مش و خطوط کانتور رنگی ارائه می شود.
فایل های ورودی
شما می توانید فایل های ورودی این آموزش را در این لینک دانلود کنید.
کد پایتون:
واردات کتابخانه های مورد نیاز
این آموزش نیاز به برخی از کتابخانه های هسته پایتون، Numpy، Matplotlib و Flopy دارد.
import os
import numpy as np
import matplotlib.pyplot as plt
from collections import OrderedDict
import flopy
from flopy.utils.triangle import Triangle as Triangle
flopy is installed in E:\Software\Anaconda3\lib\site-packages\flopy
مسیر پرونده مدل و اجرایی را تعریف کنید
پوشه برای فایل های مدل و خروجی شبیه سازی نیز تعریف شده است و همچنین اجرایی برای ژنراتور مش مثلثی و Modflow 6 تعریف شده است.
workspace = '../Model/'
triExeName = '../Exe/triangle.exe'
mf6ExeName = '../Exe/mf6.exe'
تعریف دامنه و محل تخلیه
دامنه فعال مدل یک trapezoid با سمت موازی کوچکتر است که در مدل شرقی قرار دارد. تخلیه به عنوان یک مستطیل در مرکز سمت چپ مدل مورد مفهوم قرار گرفت.
active_domain = [(0, 0), (100, 10), (100, 55), (0, 65)]
drain_domain= [(80, 31), (95, 31), (95, 34), (80, 34)]
ساختن Tringle mesh
ما از ابزار Triangle از utils flopy استفاده می کنیم که در واقع برنامه Triangle برای تولید مش ها را اجرا می کند. این آموزش همراه با اجرایی مثلث برای ویندوز همراه است، اگر شما از یک سیستم عامل دیگر استفاده می کنید که لازم است کامپایل کد منبع را که می توانید از آن پیدا کنید
https://www.cs.cmu.edu/~quake/triangle.html.
tri = Triangle(angle=30, model_ws=workspace, exe_name=triExeName)
tri.add_polygon(active_domain)
tri.add_polygon(drain_domain)
tri.add_region((10,10),0,maximum_area=30) #coarse discretization
tri.add_region((88,33),1,maximum_area=5) #fine discretization
tri.build()
تولید مثلثی و نمایندگی مرزی
گام بعدی بازنویسی مش و شناسایی مرزهای ایجاد شده از نسل مش است.
fig = plt.figure(figsize=(15, 15))
ax = plt.subplot(1, 1, 1, aspect='equal')
tri.plot(edgecolor='gray')
<matplotlib.collections.PatchCollection at 0x9005c88>
### Indentification and representation of grid boundaries
fig = plt.figure(figsize=(15, 15))
ax = plt.subplot(1, 1, 1, aspect='equal')
tri.plot(edgecolor='gray')
numberBoundaries = tri.get_boundary_marker_array().max()+1
cmap = plt.cm.get_cmap("hsv", numberBoundaries)
labelList = list(range(1,numberBoundaries))
i = 0
for ibm in range(1,numberBoundaries):
tri.plot_boundary(ibm=ibm, ax=ax,marker='o', color=cmap(ibm), label= ibm)
handles, labels = plt.gca().get_legend_handles_labels()
by_label = OrderedDict(zip(labels, handles))
plt.legend(by_label.values(), by_label.keys(), loc='lower right')
<matplotlib.legend.Legend at 0x91c6ac8>
نمایندگی ویژگی های مثلث
با استفاده از روش های موجود برای شیء مثلث "tri" ما می توانیم رأس و centroid را با موقعیت و برچسب خود نشان دهیم.
fig = plt.figure(figsize=(15, 15))
ax = plt.subplot(1, 1, 1, aspect='equal')
tri.plot(ax=ax, edgecolor='gray')
tri.plot_vertices(ax=ax, marker='o', color='blue')
tri.label_vertices(ax=ax, fontsize=10, color='blue')
tri.plot_centroids(ax=ax, marker='o', color='red')
tri.label_cells(ax=ax, fontsize=10, color='red')
output_12_0.png
پیکربندی مدل MODFLOW 6
هنگامی که مش با مثلثی داریم، می توانیم مدل MODFLOW 6 را ایجاد کنیم. همانطور که می دانید، MODFLOW 6 بین مدل و شبیه سازی تفاوت دارد، زیرا مفهوم مبادله بین مدل ها را پیاده سازی می کند. یک شبیه سازی می تواند مدل های زیادی داشته باشد، برای این مورد ما یک شبیه سازی فقط یک مدل را داریم.
name = 'mf'
sim = flopy.mf6.MFSimulation(sim_name=name, version='mf6',
exe_name=mf6ExeName,
sim_ws=workspace)
tdis = flopy.mf6.ModflowTdis(sim, time_units='SECONDS',
perioddata=[[1.0, 1, 1.]])
gwf = flopy.mf6.ModflowGwf(sim, modelname=name, save_flows=True)
ims = flopy.mf6.ModflowIms(sim, print_option='SUMMARY', complexity='complex',
outer_hclose=1.e-8, inner_hclose=1.e-8)
پیکربندی بسته بندی اختیاری مثلثی
MODFLOW 6 فایل فشرده سازی فضایی و فایل فشرده سازی زمانی را متفاوت می کند. در گزینه های اختیاری فضایی ما سه نوع شبکه داریم: منظم (.dis)، مثلثی (.disv) و بدون ساختار (* .disu). این مورد مورد مطالعه از سه جهته استفاده می کند.
cell2d = tri.get_cell2d()
vertices = tri.get_vertices()
xcyc = tri.get_xcyc()
nlay = 1
ncpl = tri.ncpl
nvert = tri.nvert
top = 1.
botm = [0.]
dis = flopy.mf6.ModflowGwfdisv(gwf, length_units='METERS',
nlay=nlay, ncpl=ncpl, nvert=nvert,
top=top, botm=botm,
vertices=vertices, cell2d=cell2d)
Configuration of the flow package and the initial conditions
npf = flopy.mf6.ModflowGwfnpf(gwf, k=0.0001)
ic = flopy.mf6.ModflowGwfic(gwf, strt=10.0)
تعریف مرز ثابت CHD مرز شرایط
طرف چپ و راست مدل به عنوان مرز CHD تعریف شده است که نشان دهنده یک جریان منطقه ای به راست از سمت هیدرولیک 30 تا 10 متر است. ببینید که لبه های لبه از لبه های شبکه مدل بدست می آید.
chdList = []
leftcells = tri.get_edge_cells(4)
rightcells = tri.get_edge_cells(2)
for icpl in leftcells: chdList.append([(0, icpl), 30])
for icpl in rightcells: chdList.append([(0, icpl), 10])
chd = flopy.mf6.ModflowGwfchd(gwf, stress_period_data=chdList)
تعریف شرایط مرزی تخلیه DRN
در قسمت سمت چپ بخش مدل، یک تخلیه به عنوان دو خط در مرز از یک ناحیه اعلام شده وجود دارد. تخلیه ارتفاع 2.5 متر پایین تر از مرز ثابت سر شرقی است.
drnList = []
drnCells = tri.get_edge_cells(5) + tri.get_edge_cells(7)
for icpl in drnCells: drnList.append([(0, icpl), 7.5, 0.10])
chd = flopy.mf6.ModflowGwfdrn(gwf, stress_period_data=drnList)
گزینه های خروجی را تعریف کنید، فایل های ورودی را وارد کنید و مدل را اجرا کنید
این آموزش نوع داده ای که ذخیره می شود و سوابق در پرونده لیست تعیین می کند. سپس فایل های مدل Modflow 6 روی پوشه مدل نوشته شده و شبیه سازی انجام می شود.
oc = flopy.mf6.ModflowGwfoc(gwf,
budget_filerecord='{}.cbc'.format(name),
head_filerecord='{}.hds'.format(name),
saverecord=[('HEAD', 'LAST'),
('BUDGET', 'LAST')],
printrecord=[('HEAD', 'LAST'),
('BUDGET', 'LAST')])
sim.write_simulation()
success, buff = sim.run_simulation()
writing simulation...
writing simulation name file...
writing simulation tdis package...
writing ims package ims_-1...
writing model mf...
writing model name file...
writing package disv...
writing package npf...
writing package ic...
writing package chd_0...
INFORMATION: maxbound in ('gwf6', 'chd', 'dimensions') changed to 19 based on size of stress_period_data
writing package drn_0...
INFORMATION: maxbound in ('gwf6', 'drn', 'dimensions') changed to 32 based on size of stress_period_data
writing package oc...
FloPy is using the following executable to run the model: ../Exe/mf6.exe
MODFLOW 6
U.S. GEOLOGICAL SURVEY MODULAR HYDROLOGIC MODEL
VERSION 6.0.3 08/09/2018
MODFLOW 6 compiled Aug 09 2018 13:40:32 with IFORT compiler (ver. 18.0.3)
This software has been approved for release by the U.S. Geological
Survey (USGS). Although the software has been subjected to rigorous
review, the USGS reserves the right to update the software as needed
pursuant to further analysis and review. No warranty, expressed or
implied, is made by the USGS or the U.S. Government as to the
functionality of the software and related material nor shall the
fact of release constitute any such warranty. Furthermore, the
software is released on condition that neither the USGS nor the U.S.
Government shall be held liable for any damages resulting from its
authorized or unauthorized use. Also refer to the USGS Water
Resources Software User Rights Notice for complete use, copyright,
and distribution information.
Run start date and time (yyyy/mm/dd hh:mm:ss): 2019/03/27 10:24:15
Writing simulation list file: mfsim.lst
Using Simulation name file: mfsim.nam
Solving: Stress period: 1 Time step: 1
Run end date and time (yyyy/mm/dd hh:mm:ss): 2019/03/27 10:24:15
Elapsed run time: 0.041 Seconds
خاتمه طبیعی شبیه سازی.
نمایندگی داده های خروجی
ما تراز محاسبه شده را وارد می کنیم و آنها را در فرمت مدل نشان می دهیم. نمایندگی از Equipotential تراز با Interpolation تراز در هر centocroid سلول انجام می شود.
fname = os.path.join(workspace, name + '.hds')
hdobj = flopy.utils.HeadFile(fname, precision='double')
head = hdobj.get_data()
fig = plt.figure(figsize=(15, 15))
ax = plt.subplot(1, 1, 1, aspect='equal')
img=tri.plot(ax=ax, a=head[0, 0, :], cmap='Spectral')
fig.colorbar(img, fraction=0.02)
<matplotlib.colorbar.Colorbar at 0xa6eacf8>
from scipy.interpolate import griddata
x, y = tri.get_xcyc()[:,0],tri.get_xcyc()[:,1]
xG = np.linspace(x.min(),x.max(),100)
yG = np.linspace(y.min(),y.max(),100)
X, Y = np.meshgrid(xG, yG)
z = head[0,0,:]
fig = plt.figure(figsize=(15, 15))
ax = plt.subplot(1, 1, 1, aspect='equal')
img=tri.plot(ax=ax, edgecolor='lightsteelblue', a=head[0, 0, :], cmap='Spectral', alpha=0.8)
fig.colorbar(img, fraction=0.02)
Ti = griddata((x, y), z, (X, Y), method='cubic')
CS = ax.contour(X,Y,Ti, colors='Cyan', linewidths=2.0)
ax.clabel(CS, inline=1, fontsize=15, fmt='%1.0f')
plt.show()
برای یافتن تمامی مطالب مرتبط با این مطلب در سایت از جستجوی سایت در حاشیه سمت راست و بالای صفجه استفاده فرمایید.
ورود به بخش آموزش های متنی GMS
دانلود آخرین نسخه نرم افزار GMS
دریافت لایسنس ارزیابی (14 روزه)
برای سفارش انجام مدل سازی اینجا کلیک کنید
شناسه تلگرام مدیر سایت: SubBasin@
نشانی ایمیل: behzadsarhadi@gmail.com
(سوالات تخصصی را در گروه تلگرام ارسال کنید)
_______________________________________________________
نظرات (۰)