とらりもん - 開水路の水面形 Diff
- Added parts are displayed like this.
- Deleted parts are displayed
like this.
!!問題設定
以下の図のように、水路を水が横方向に流れている状況を考えよう。川幅(水路幅)は一定とするので、奥行き方向は考えない。上流からの給水と下流での排水が常に一定で等しい状態、つまり流れが定常状態にあるとする。
{{attach_view(concept.png)}}
こういう状況で、任意に河床の形状z(x)が与えられたときに、水面形がどうなるか、計算しよう。
!!基礎方程式
この問題を扱う出発点として, 以下のような2つの基礎方程式がある:
!!!水面近傍におけるエネルギー保存則(ベルヌーイの定理):
{{attach_view(eq1.png)}}
.............. (1)
ここで、Uはエネルギー(一定), ρは水の密度、gは重力加速度、hは水面の高さ、vは流速
!!!質量保存則(流量一定):
{{attach_view(eq2.png)}}
..................(2)
ここで、qは単位幅流量(一定)、zは河床の高さ
以上の2つの式が、基礎方程式である。
!!微分方程式
(1)と(2)から、以下の式が導かれる(水面と流速に関する連立微分方程式):
{{attach_view(eq3.png)}}
......(3)
これが、解くべき微分方程式である。
ちなみに、フルード数Frを導入して(3)を変形すると、
{{attach_view(eq4.png)}}
............(4)
のようになる。これをみると、フルード数が1に近付くと右辺が発散することがわかる。従って、フルード数が1をまたいで変化するような場合(跳水など)はこの方程式は解けない。
問: (1)(2)より、(3)を導け。
ヒント: Uやqは、場所xによらない定数なので、xで微分すればゼロになって消える。
問: (3)より、(4)を導け。
数値解析
式(3)は、1階の連立常微分方程式なので、これまで学んだ方法で数値解析できる。
問: 式(3)を、以下の条件で、ルンゲ・クッタ法で解き、結果をグラフにせよ。
h(0)=0.5 [m]
v(0)=1.0 [m/sec]
z(x)=0.1sin(x) [m]
ρ=1000 [kg/m3]
g=9.8 [m/s2]
xは0から10 [m]
ヒント: 河床形状を表現するために、三角関数を使う必要がある。三角関数のような解析関数をC言語で使うには、数学ライブラリを利用する。具体的には、# include <stdio.h>の次の行に、# include <math.h>という行を入れる。コンパイルするときは、gccに-lmというオプションをつけること。
プログラムの例
出力例:
{{attach_view(case1.png)}}
参考: これを描いたGnuplotコマンド
set xlabel [m]
set ylabel [m]
set yrange [-0.3:1]
plot "out" using 1:2 title "bed" w l, \
"out" using 1:3 title "water surface" w l 3
問: この解のフルード数もプロットしてみよ。
!!常流と射流
上の条件(常流)と同じエネルギーおよび同じ流量を実現するような、別の解(射流)も存在するはずである。それを見つけてみよう。
問: まず、水深と流速に対してエネルギーを与える式、つまり式(1)を変形して、水深と流量に対してエネルギーを与える式を作れ。ヒント:式(2)を用いて、vをqとhで表す。平均河床を考えて、z=0としてよい。
問: その式を、横軸を水深として縦軸をエネルギーとするグラフにプロットしてみよ。その際、密度ρを入れたままにすると、SI単位系ではエネルギーUが莫大な数字になるので、エネルギーUのかわりに、U/ρを縦軸にとるとよい。
{{attach_view(energy.png)}}
参考: これを描いたgnuplotコマンド
set xrange [0:1.0]
set yrange [0:12]
set xlabel "water depth [m]"
set ylabel "energy [m^2/s^2]"
set xtic 0.1
set ytic 1
plot 9.8*x title "gh", 0.5*0.5/(x*x*2) title "q^2/(2h^2)", 9.8*x+0.5*0.5/(x*x*2) title "U"
これを見ると、水深0.5 mに対応するエネルギー(縦軸で5.5のあたり)は、水深約0.18 mに対応するエネルギーとおなじくらいであることがわかる。つまり、上の条件で計算したのと同じエネルギーかつ同じ流量を、水深0.18 m程度の流れでも実現できることがわかる。(これが射流)
{{attach_view(energy2.png)}}
この射流に対応する(左端での)水深を厳密に求めるには、上で導いたUの式に関する代数方程式を解けば良い。それには、Newton法などを用いる。
問: この射流に対応する(左端での)水深をNewton法で求めよ。
プログラムの例:
/* openchannel.c
* 2004/01/26 Kenlo Nishida
* compile: $ gcc openchannel.c -o openchannel -lm
* usage: Give h and v from stdin.
*/
# include <stdio.h>
# include <stdlib.h>
# include <math.h>
# define NV 2 // number of (dependent) variables for Runge-Kutta solver
# define DX 0.01 // step along x-direction
# define G 9.8 // gravity acceleration
# define rho 1000 // water density [kg/m^3]
double z(double x){ // river bed
return(0.1*sin(x));}
double dzdx(double x){ // gradient (slope) of the river bed, i.e., dz/dx
return ( (z(x+DX)-z(x-DX))/(2*DX) );}
void func(double x, double *p, double *dpdx) // differential equation
// x: x-coordinate
// p[0] ... h
// p[1] ... v
// dpdx is output
{double h, v;
h=p[0]; v=p[1];
dpdx[0] = dzdx(x)*v*v/(v*v - (h-z(x))*G); // dh/dx
dpdx[1] = -1*dzdx(x)*v*G/(v*v - (h-z(x))*G); // dv/dx
}
void RungeKutta(double x, double dx, double *p){
double dpdx[NV], dpdx0[NV], dpdx1[NV],dpdx2[NV],dpdx3[NV];
double p1[NV], p2[NV], p3[NV];
int i;
func(x, p, dpdx0); for (i=0; i<NV; i++) p1[i]=p[i]+dpdx0[i]*dx/2;
func(x+dx/2, p1, dpdx1); for (i=0; i<NV; i++) p2[i]=p[i]+dpdx1[i]*dx/2;
func(x+dx/2, p2, dpdx2); for (i=0; i<NV; i++) p3[i]=p[i]+dpdx2[i]*dx;
func(x+dx, p3, dpdx3);
for (i=0; i<NV; i++) dpdx[i]=(dpdx0[i] + 2*dpdx1[i] + 2*dpdx2[i] + dpdx3[i])/6;
for (i=0; i<NV; i++) p[i]=p[i]+dpdx[i]*dx;
}
int main(){
double x;
double x0=0, x1=10; // boundary
double p[NV];
double dx=DX;
// setting initial condition
scanf("%lf", &p[0]); // h
scanf("%lf", &p[1]); // v
for (x=x0; x<x1; x=x+dx){
printf("%lf %lf %lf %lf\n", x, z(x), p[0], p[1]/sqrt(G*(p[0]-z(x))));
RungeKutta(x, dx, p);
}
}
その結果、h = 0.187 mが、この射流の条件であることがわかる。
問: この射流(左端での水深は0.187 m)を、上と同様に計算し、プロットしてみよ。
参考になるページ: [[開水路 (Wikipedia)|https://ja.wikipedia.org/wiki/%E9%96%8B%E6%B0%B4%E8%B7%AF]]
以下の図のように、水路を水が横方向に流れている状況を考えよう。川幅(水路幅)は一定とするので、奥行き方向は考えない。上流からの給水と下流での排水が常に一定で等しい状態、つまり流れが定常状態にあるとする。
{{attach_view(concept.png)}}
こういう状況で、任意に河床の形状z(x)が与えられたときに、水面形がどうなるか、計算しよう。
!!基礎方程式
この問題を扱う出発点として, 以下のような2つの基礎方程式がある:
!!!水面近傍におけるエネルギー保存則(ベルヌーイの定理):
{{attach_view(eq1.png)}}
.............. (1)
ここで、Uはエネルギー(一定), ρは水の密度、gは重力加速度、hは水面の高さ、vは流速
!!!質量保存則(流量一定):
{{attach_view(eq2.png)}}
..................(2)
ここで、qは単位幅流量(一定)、zは河床の高さ
以上の2つの式が、基礎方程式である。
!!微分方程式
(1)と(2)から、以下の式が導かれる(水面と流速に関する連立微分方程式):
{{attach_view(eq3.png)}}
......(3)
これが、解くべき微分方程式である。
ちなみに、フルード数Frを導入して(3)を変形すると、
{{attach_view(eq4.png)}}
............(4)
のようになる。これをみると、フルード数が1に近付くと右辺が発散することがわかる。従って、フルード数が1をまたいで変化するような場合(跳水など)はこの方程式は解けない。
問: (1)(2)より、(3)を導け。
ヒント: Uやqは、場所xによらない定数なので、xで微分すればゼロになって消える。
問: (3)より、(4)を導け。
数値解析
式(3)は、1階の連立常微分方程式なので、これまで学んだ方法で数値解析できる。
問: 式(3)を、以下の条件で、ルンゲ・クッタ法で解き、結果をグラフにせよ。
h(0)=0.5 [m]
v(0)=1.0 [m/sec]
z(x)=0.1sin(x) [m]
ρ=1000 [kg/m3]
g=9.8 [m/s2]
xは0から10 [m]
ヒント: 河床形状を表現するために、三角関数を使う必要がある。三角関数のような解析関数をC言語で使うには、数学ライブラリを利用する。具体的には、# include <stdio.h>の次の行に、# include <math.h>という行を入れる。コンパイルするときは、gccに-lmというオプションをつけること。
プログラムの例
出力例:
{{attach_view(case1.png)}}
参考: これを描いたGnuplotコマンド
set xlabel [m]
set ylabel [m]
set yrange [-0.3:1]
plot "out" using 1:2 title "bed" w l, \
"out" using 1:3 title "water surface" w l 3
問: この解のフルード数もプロットしてみよ。
!!常流と射流
上の条件(常流)と同じエネルギーおよび同じ流量を実現するような、別の解(射流)も存在するはずである。それを見つけてみよう。
問: まず、水深と流速に対してエネルギーを与える式、つまり式(1)を変形して、水深と流量に対してエネルギーを与える式を作れ。ヒント:式(2)を用いて、vをqとhで表す。平均河床を考えて、z=0としてよい。
問: その式を、横軸を水深として縦軸をエネルギーとするグラフにプロットしてみよ。その際、密度ρを入れたままにすると、SI単位系ではエネルギーUが莫大な数字になるので、エネルギーUのかわりに、U/ρを縦軸にとるとよい。
{{attach_view(energy.png)}}
参考: これを描いたgnuplotコマンド
set xrange [0:1.0]
set yrange [0:12]
set xlabel "water depth [m]"
set ylabel "energy [m^2/s^2]"
set xtic 0.1
set ytic 1
plot 9.8*x title "gh", 0.5*0.5/(x*x*2) title "q^2/(2h^2)", 9.8*x+0.5*0.5/(x*x*2) title "U"
これを見ると、水深0.5 mに対応するエネルギー(縦軸で5.5のあたり)は、水深約0.18 mに対応するエネルギーとおなじくらいであることがわかる。つまり、上の条件で計算したのと同じエネルギーかつ同じ流量を、水深0.18 m程度の流れでも実現できることがわかる。(これが射流)
{{attach_view(energy2.png)}}
この射流に対応する(左端での)水深を厳密に求めるには、上で導いたUの式に関する代数方程式を解けば良い。それには、Newton法などを用いる。
問: この射流に対応する(左端での)水深をNewton法で求めよ。
プログラムの例:
/* openchannel.c
* 2004/01/26 Kenlo Nishida
* compile: $ gcc openchannel.c -o openchannel -lm
* usage: Give h and v from stdin.
*/
# include <stdio.h>
# include <stdlib.h>
# include <math.h>
# define NV 2 // number of (dependent) variables for Runge-Kutta solver
# define DX 0.01 // step along x-direction
# define G 9.8 // gravity acceleration
# define rho 1000 // water density [kg/m^3]
double z(double x){ // river bed
return(0.1*sin(x));}
double dzdx(double x){ // gradient (slope) of the river bed, i.e., dz/dx
return ( (z(x+DX)-z(x-DX))/(2*DX) );}
void func(double x, double *p, double *dpdx) // differential equation
// x: x-coordinate
// p[0] ... h
// p[1] ... v
// dpdx is output
{double h, v;
h=p[0]; v=p[1];
dpdx[0] = dzdx(x)*v*v/(v*v - (h-z(x))*G); // dh/dx
dpdx[1] = -1*dzdx(x)*v*G/(v*v - (h-z(x))*G); // dv/dx
}
void RungeKutta(double x, double dx, double *p){
double dpdx[NV], dpdx0[NV], dpdx1[NV],dpdx2[NV],dpdx3[NV];
double p1[NV], p2[NV], p3[NV];
int i;
func(x, p, dpdx0); for (i=0; i<NV; i++) p1[i]=p[i]+dpdx0[i]*dx/2;
func(x+dx/2, p1, dpdx1); for (i=0; i<NV; i++) p2[i]=p[i]+dpdx1[i]*dx/2;
func(x+dx/2, p2, dpdx2); for (i=0; i<NV; i++) p3[i]=p[i]+dpdx2[i]*dx;
func(x+dx, p3, dpdx3);
for (i=0; i<NV; i++) dpdx[i]=(dpdx0[i] + 2*dpdx1[i] + 2*dpdx2[i] + dpdx3[i])/6;
for (i=0; i<NV; i++) p[i]=p[i]+dpdx[i]*dx;
}
int main(){
double x;
double x0=0, x1=10; // boundary
double p[NV];
double dx=DX;
// setting initial condition
scanf("%lf", &p[0]); // h
scanf("%lf", &p[1]); // v
for (x=x0; x<x1; x=x+dx){
printf("%lf %lf %lf %lf\n", x, z(x), p[0], p[1]/sqrt(G*(p[0]-z(x))));
RungeKutta(x, dx, p);
}
}
その結果、h = 0.187 mが、この射流の条件であることがわかる。
問: この射流(左端での水深は0.187 m)を、上と同様に計算し、プロットしてみよ。
参考になるページ: [[開水路 (Wikipedia)|https://ja.wikipedia.org/wiki/%E9%96%8B%E6%B0%B4%E8%B7%AF]]