GMS建模
===
建立XYZ檔
---
**1. 開QGIS將所選DEM圖層轉成.XYZ檔**
.xyz檔就是將DTM每個點位轉成X軸、Y軸、Z軸顯示
內政部20M DTM 使用的CRS(大地地理座標) 為 TWD97-121 (EPSG:3826)
**2. .XYZ檔打開後看到第一行(286XXX)為經度(X軸)
第二行(276XXX)為緯度(Y軸)
第三行為高程(Z軸)**
**3. 將.XYZ檔用Excel打開編輯，將X、Y軸從原點(0,0)開始重新編輯**
**4. 排序X軸順序
設置表格後 > 檢視 > 凍結表格 > 凍結首欄
常用 > 排序與篩選 > 由小至大**
**5. 記住有多少筆資料**
如果有一行是XYZ軸記得去掉
**6. 將排序完的資料另存新檔(.xyz)**
**7. 最後加上一些東西即可**
```
TIN
BEGT
ACTIVETIN
ID 1
TNAM "20mtest" (檔案名稱)
THORID "0"
MAT 1
GUID 978bf581-e7f2-4894-954e-eb7ba01da9b0 (識別碼)
VERT 18733(第五點中的資料筆數)
X Y Z數據...
TRI 0
ENDT
```
**P.S. 20m DEM XY以20為等距 6m DEM 以6為等距**
輸入檔案
---
**開啟GMS > 左上File > Open 打開.xyz檔 > 右鍵Tin Data裡的檔案 > Z高程建立資料庫**
**從上方工具列找Display Options > TIN Data > 勾選Vertices顯示DTM節點**
**TINs > Triangulate將各節點相連形成不規則三角形 (點→線)**
建立2D Mesh Data
---
**在上方 materals建立分層資料還有顏色(ex.崩積層,泥岩,砂岩...等)**
**TINs > Convert to > 2D Mesh 將不規則三角形節點轉成網格 (線→面)**
**最上層地層顏色為網格顏色**
設置鑽井資料
---
**左邊工作欄右鍵 New > Borehole 建立鑽井資料**
**要移至平面才能開始操作**
**建設所要得井位**
**編輯井位資料**
根據設置井位的XY軸，透過之前的XYZ檔的高程去假設出地層
**Display Option > Borehole Data 去調整井位顯示**
**全部井位確定建好 > 上方Borehole > Auto-Assign Horizons將各地層作編號分類**
建立3D Mesh Data
---
**至少要有==6個==地下水井且分散面的四周(不能太集中)才能建立3D Mesh**
**頂層使用Tin data，Constant自訂義模型底部高程多少**
**生成3維網格細節**
* Interpolation Method (使用的算法)
1. Inverse Distance Weighted (反距離權重法)
2. Natural Neighbor (自然對數法)
真實井位多用反距離，真實井位少用自然對數
* Nodel Function (節點內插法)
1. Constant (根據各點位置內插)
2. Gradient Plane (根據頂層坡度內插)
3. Quadratic (根據網格密度品質內插)
表面水平用 Constant ， 表面起伏用 Gradient
* Meshing Options (網格型式)
1. Prevent Pinchouts (只生成Pyramids & Hexahedra)
2. Minimum Element Thickness (網格最小高度)
3. Refine Elements (固定網格型態)
a. Refine All Elements (全部模型只使用 四面體網格Tetrahedral)
b. Refine Layers With Pyramid Elements (全部模型只使用 五面體網格Pyramid)
* Subdivide Material Layers (各地層分層)
1. Target Layer Thickness (每個網格固定高度)
2. Number of Layers (每個地層分幾層)
**跑出3D Mesh圖**
單純看建模情形可分各地層1層來看，要導出確定模型再細部分多層
不然一次分太多層對電腦記憶體造成大量負擔
模型處理&輸出
---
**右鍵3D Mesh裡檔案 > 輸出3D Mesh (mesh1.txt)**
**用Excel打開剛導出的檔案 (mesh1.txt)**
Excel有資料筆數限制，如果資料超過一百萬筆，以每一百萬為界線將資料分割後處理。
**1. 將E4T與ND分開個別處理**
**2. E4T刪除,存e4t (.csv)**
**3. ND刪除,存node(.csv)**
**4. 將<font color="#8F16C6">GMS2Flac3D</font> (Matlab Code)放進e4t(.csv)與node(.csv)同個資料夾並執行**
==**以下為全四面體轉換程式碼**== (有勾選Refine element)
```
function GMS2Flac3D
%從flac3d中取得資訊
tic
[xyz,elem]=read_GMS_files;
disp('GMS files read...')
toc
[k]=read(elem,xyz);
toc
save workspace
%diary('D:\FORMAT\lol.txt');diary on;
%% *******************************
%% * FUNCTIONS *
%% *******************************
function [xyz,elem]=read_GMS_files
fprintf('node.csv：\n');
type node.csv
xyz = csvread('node.csv');
fprintf('e4t.csv：\n');
type e4t.csv
elem = csvread('e4t.csv');
end
function [k]=read(elem,xyz)
fid=fopen('GMS2FLAC.flac3d','a+');
fprintf(fid,'* FLAC3D grid produced by FLAC3D\n');
fprintf(fid,'* GRIDPOINTS\n');
[m,~]=size(xyz);
for i=1:m
nn=xyz(i,1);
xx=xyz(i,2);
yy=xyz(i,3);
zz=xyz(i,4);
fprintf(fid,'G %d %e %e %e\n',nn,xx,yy,zz);
end
fprintf(fid,'* ZONES\n');
[n,~]=size(elem);
for i=1:n
fN=elem(i,1);
f0=elem(i,2);
f1=elem(i,3);
f2=elem(i,4);
f3=elem(i,5);
fprintf(fid,'Z T4 %d %d %d %d %d\n',fN,f0,f1,f2,f3);
clear n
end
fprintf(fid,'* GROUPS\n');
qq=6; %NUMBER OF GROUPS
[h,~]=size(elem);
for i=1:qq
G=i;
fprintf(fid,'ZGROUP %d\n',G);
for j=1:h
k=[];
GR=num2str(elem(j,6));
bb=num2str(i);
if GR==bb
w=(elem(j,1));
k=(w);
fprintf(fid,'%d\n',k);
j=j+1;
else
disp('false')
end
end
end
clear h
end
end
```
==**以下為四面、五面、六面體轉換程式**== (無勾選Refine element)
五面體與六面體網格處理方法同E4T
```
function GMS2Flac3D
%從flac3d中取得資訊
tic
[ND,e4t,e5p,e6w]=read_GMS_files;
disp('GMS files read...')
toc
[k]=read(e4t,e5p,e6w,ND);
toc
save workspace
%% *******************************
%% * FUNCTIONS *
%% *******************************
function [ND,e4t,e5p,e6w]=read_GMS_files
fprintf('node.csv：\n');
type node.csv
ND = csvread('node.csv');
fprintf('e4t.csv：\n');
type e4t.csv
e4t = csvread('e4t.csv');
fprintf('e5p.csv：\n');
type e5p.csv
e5p = csvread('e5p.csv');
fprintf('e6w.csv：\n');
type e6w.csv
e6w = csvread('e6w.csv');
end
function [k]=read(e4t,e5p,e6w,ND)
fid=fopen('GMS2FLAC.flac3d','a+');
fprintf(fid,'* FLAC3D grid produced by FLAC3D\n');
fprintf(fid,'* GRIDPOINTS\n');
[gridpoint,~]=size(ND);
for i=1:gridpoint
nn=ND(i,1);
xx=ND(i,2);
yy=ND(i,3);
zz=ND(i,4);
fprintf(fid,'G %d %e %e %e\n',nn,xx,yy,zz);
clear gridpoint
end
fprintf(fid,'* ZONES\n');
[tetrahedral,~]=size(e4t);
for i=1:tetrahedral
fN=e4t(i,1); %ID
f0=e4t(i,2);
f1=e4t(i,3);
f2=e4t(i,4);
f3=e4t(i,5);
fprintf(fid,'Z T4 %d %d %d %d %d\n',fN,f0,f1,f2,f3);
clear tetrahedral
end
[pyramid,~]=size(e5p);
for i=1:pyramid
fN=e5p(i,1); %ID
f0=e5p(i,2);
f1=e5p(i,3);
f2=e5p(i,4);
f3=e5p(i,5);
f4=e5p(i,6);
fprintf(fid,'Z P5 %d %d %d %d %d %d\n',fN,f0,f1,f3,f4,f2);
clear pyramid
end
[wedge,~]=size(e6w);
for i=1:wedge
fN=e6w(i,1); %ID
f0=e6w(i,2);
f1=e6w(i,3);
f2=e6w(i,4);
f3=e6w(i,5);
f4=e6w(i,6);
f5=e6w(i,7);
fprintf(fid,'Z W6 %d %d %d %d %d %d %d\n',fN,f0,f2,f3,f1,f5,f4);
clear wedge
end
fprintf(fid,'* GROUPS\n');
strata=4; %地層個數
[t4f,~]=size(e4t);
[p5e,~]=size(e5p);
[w6e,~]=size(e6w);
for i=1:strata
G=i;
fprintf(fid,'ZGROUP %d\n',G);
for j=1:t4f
k=[];
GR=num2str(e4t(j,6));
bb=num2str(i);
if GR==bb
w=(e4t(j,1));
k=(w);
fprintf(fid,'%d\n',k);
j=j+1;
else
disp('false')
end
end
for j=1:p5e
k=[];
GR=num2str(e5p(j,7));
bb=num2str(i);
if GR==bb
w=(e5p(j,1));
k=(w);
fprintf(fid,'%d\n',k);
j=j+1;
else
disp('false')
end
end
for j=1:w6e
k=[];
GR=num2str(e6w(j,8));
bb=num2str(i);
if GR==bb
w=(e6w(j,1));
k=(w);
fprintf(fid,'%d\n',k);
j=j+1;
else
disp('false')
end
end
end
clear t4f
clear p5e
clear w6e
end
end
```
**5. 跑完會生成<font color="#8F16C6">GMS2FLAC.flac3D</font>(如果要重跑要刪除後再跑)**
模型導入Flac3D
---
**將生成的GMS2FLAC.flac3d 用Flac3D 打開便能看見建立的模型**
左上File > Open > GMS2FLAC.flac3d
**打開Model視窗可以看到自己建立的模型**
**點選地層 > 工具列 Assign a Group... > 將所有地層做分類**
**工具列Assign Group names to face... > Ignore existing group names分東西南北面**
分面是要方便後面Flac3d做力學運算
**回到Plot視窗 > 右上+Zone可以看見分完地層的模型**
**==看到模型沒問題在下方輸入 model save 'ini_model' 將模型存檔==**
生成地下水位面
---
**1. 在borehole裡輸入地下水位高程**
**2. borehole(上面) > 選advanced > waterable... (幫他取名字water之類的)**
**3. 對Scatter右鍵 > interpolate to > active tin (轉成.tin檔,所以要在Tin Data裡面看)**
**4. 將地下水位.tin檔 > 生成2D Mesh**
**5. 2D Mesh > 生成 Cad圖檔(.dxf) > 將Cad圖檔存至運算的Flac3D資料夾內**
**6. 透過Flac3D程式碼便可在Water Table視窗看到地下水位面**
地下水位面(舊方法)
---
**4. 接續上面第3點到Tin Data裡將地層與地下水位資料輸出**
要輸出地層先點default再點surface輸出.tin檔
要輸出水位先點watertable再點surface輸出.tin檔
將輸出的地下水位面資料(surface.tin)用Matlab轉換成Flac3D可讀取的格式
```
clc
clear
X0_=0; %X軸原點
Y0_=0; %Y軸原點
BOT=50; %Z軸原點
xsize_=615; %X軸邊長
ysize_=355; %Y軸邊長
interval_=5; %節點間距(網格大小,DTM精度)
yline_=(ysize_/interval_)+1;
xline_=(xsize_/interval_)+1;
todata=yline_*xline_;
load('surface.tin');
XX=surface(:,1);
YY=surface(:,2);
TT=surface(:,3);
X1=X0_:interval_:xsize_;
Y1=Y0_:interval_:ysize_;
[X2,Y2]=meshgrid(X1,Y1);
Z1=griddata(XX,YY,TT,X2,Y2);
surf(X1,Y1,Z1)
%將文件轉換成Flac3D
for i=1:yline_
no_(i,1:xline_)=i+100;
BOTNO(i,1:xline_)=i+100;
BOT_(i,1:xline_)=BOT;
end
DATA_(:,1)=flipud(rot90(reshape(rot90(no_), 1, todata)));
XD=flipud(rot90(X2));
ZD=flipud(rot90(Z1));
DATA_(:,2)=reshape(XD, 1, todata);
DATA_(:,3)=reshape(ZD, 1, todata);
DATA_(isnan(DATA_(:,3)),3)=-0.5;
DATAB_(:,1)=flipud(rot90(reshape(rot90(BOTNO), 1, todata)));
DATAB_(:,2)=DATA_(:,2);
DATAB_(:,3)=reshape(BOT_, 1, todata);
ch(1:todata,1)=char('t');
ch(1:todata,2)=char('a');
ch(1:todata,3)=char('b');
ch(1:todata,4)=char('l');
ch(1:todata,5)=char('e');
ch(1:todata,6)=char(' ');
ftable1_=num2str(DATA_(:,1));
fx1_=num2str(DATA_(:,2));
fz1_=num2str(DATA_(:,3));
ftable2_=num2str(DATAB_(:,1));
fx2_=num2str(DATAB_(:,2));
fz2_=num2str(DATAB_(:,3));
temp={ch,ch(1:todata,6),ftable1_,ch(1:todata,6),fx1_,ch(1:todata,6),fz1_};
temp=cell2mat(temp);
dlmwrite('watertab.dat',';Tables of interpolated data from contour map','delimiter','');
dlmwrite('watertab.dat',temp,'-append', 'delimiter','');
dlmwrite('watertab.dat',' ','delimiter','','-append')
```
{"metaMigratedAt":"2023-06-15T15:57:37.315Z","metaMigratedFrom":"Content","title":"GMS建模","breaks":true,"contributors":"[{\"id\":\"75864b70-695a-464c-8b75-3038be1cb3ff\",\"add\":21232,\"del\":9071}]"}