1.1蛙跳格式
蛙跳格式是一个很常用的差分格式,它在时间和空间上都采用中心差分。假设离散化如下图所示:
(2-1)
蛙跳格式是时间和空间上均为二阶精度的,三层一步显式差分格式。因而,为了计算(n+1)层的值,需要(n-1)层和n层的已知值。另外,还可以看出,后一偶数层的值实际上等于前一偶数层的值加上一个变化量,而跳过其间的奇数层,这就是蛙跳式这个名称的来历。
2.2对圣维南方程方程差分
上述格式用于圣维南方程,便得到Q和Z的显式表达式。
无侧向入流时完整偏微分方程为:
(1-2)
按式(1-1)对(1-2)进行离散可得:
(1-3)
上式就是由连续方程得到的在i点和j时间层的水位显式表达式。
从动量方程得:
(1-4)
于是,对于在i点和j时间层得流量显式表达式为:
(1-5)
阻力项也可以用下面差分形式:
这里Φ是权函数,取Φ=0.5则得:
从而得:
1.3稳定性
该格式的稳定性的必要条件:
考虑到阻力项的影响,∆t还必须满足:
其中 为特征流速, 为重力波速, 为摩阻比降。
Fortran版
!b=The width of the rectangular cross-section,unit:m
!length=Total length of the stream under consideration,unit:m
!g=Gravitational acceleration,unit:m/s^2
!deltat=Time step,unit:s
!dx=Spatial step
!rn=Roughness coefficient
!t=Total computation time
!h0=Initial water depth,unit:m
!qq0=Initial discharge
!areao=Cross Section Area, Unit: m^2
!cko=Hydromodulus
program leapfrog
use ifport
implicit none
integer::n,i,j,itert,t,start(2), finish(2), ierr
real b,a,g,deltat,dx,length,lambda,&
qq0,bb,rn,h0,Total_Time
real,dimension(:),allocatable::z0,z,z1,zh,&
q0,q,q1,qh,cko,bwo,areao
real area,ck,random1
open(1,file='dataf.txt',status='old')
open(2,file='outf.txt',status='old')
read(1,*)b,h0,length,qq0
read(1,*)g,deltat,dx,rn,t
call gettimeofday(start, ierr)
bb=length/dx
itert=ifix(bb)
n=itert
write(*,*)'n=',n
allocate(z(0:n+1),q(0:n+1),z0(0:n+1),q0(0:n+1))
!do i=1,n
! z(i)=h0
! q(i)=qq0
!end do
call random_seed()
do i=1,n
z(i)=random1(h0,h0+1.0)
q(i)=random1(qq0,qq0+1.0)
!write(*,*) z(i)
end do
allocate(cko(n),bwo(n),areao(0:n+1))
do i=1,n
cko(i)=ck(b,z(i),rn)
bwo(i)=b
areao(i)=area(b,z(i))
end do
z(0) = z(n)
z(n+1) = z(1)
q(0) = q(n)
q(n+1) = q(1)
areao(0)=areao(n)
areao(n+1)=areao(1)
do i=1,n+1
z0(i) = z(i)
q0(i) = q(i)
end do
lambda=deltat/dx
allocate(zh(n),qh(n),z1(0:n+1),q1(0:n+1))
!First step is FTBS
do i=1,n
zh(i)=z(i)-lambda*(1./bwo(i))*(q(i)-q(i-1))
end do
do i=1,n
qh(i)=(1./(g*deltat*(areao(i)*abs(q(i))/cko(i)**2)+1.))&
*(q(i)*(1.-(g*deltat*((areao(i)*abs(q(i)))/cko(i)**2)))-2.&
*lambda*(q(i)/areao(i))*(q(i)-q(i-1))+lambda&
*((q(i)**2*bwo(i)-g*areao(i)**3)/areao(i)**2)*&
(z(i)-z(i-1))+lambda*(q(i)**2/areao(i)**2)*z(i)*(areao(i)-areao(i-1)))
end do
do i=1,n
z1(i)=zh(i)
q1(i)=qh(i)
end do
z1(0) = z1(n)
z1(n+1) = z1(1)
q1(0) = q1(n)
q1(n+1) = q1(1)
do i=1,n+1
areao(i)=area(b,z1(i))
end do
do j=1,t/deltat
do i=1,n
zh(i)=z(i)-lambda*(1./bwo(i))*(q1(i+1)-q1(i-1))
end do
do i=1,n
qh(i)=(1./(g*deltat*(areao(i)*abs(q1(i))/cko(i)**2)+1.))&
*(q(i)*(1.-(g*deltat*((areao(i)*abs(q1(i)))/cko(i)**2)))-&
2.*lambda*(q1(i)/areao(i))*(q1(i+1)-q1(i-1))+lambda&
*((q1(i)**2*bwo(i)-g*areao(i)**3)/areao(i)**2)&
*(z1(i+1)-z1(i-1))+lambda*(q1(i)**2/areao(i)**2)*z1(i)*(areao(i+1)-areao(i-1)))
end do
do i=0,n+1
z(i)=z1(i)
q(i)=q1(i)
end do
do i=1,n
z1(i)=zh(i)
q1(i)=qh(i)
end do
z1(0) = z1(n)
z1(n+1) = z1(1)
q1(0) = q1(n)
q1(n+1) = q1(1)
end do
call gettimeofday(finish, ierr)
Total_Time = finish(1)-start(1)+(finish(2)-start(2))/1000000.0
write(*,*) 'Total time is',Total_time,'s'
do i=1,n
write(*,*) i,z0(i),z1(i),abs(z1(i)-z0(i)),q0(i),q1(i),abs(q1(i)-q0(i))
end do
end program
function area(b,y)
implicit none
real b,y,area
area=b*y
end function
function ck(b,y,rn)
implicit none
real ck,b,y,rn
ck=(b*y/(2*y+b))**0.667*b*y/rn
return
end function
function random1(lb,ub)
implicit none
integer i
real lb,ub
real random1
real t,len
len=ub-lb
call random_number(t)
random1=lb+len*t
return
end function
C语言版
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
double area(double b, double y)
{
double a;
a = b * y;
return a;
}
double ck(double b, double y, double rn)
{
double c,a,d;
a = b * y / (2 * y + b);
d = pow(a, 0.667);
c = d*b*y / rn;
return c;
}
void main()
{
int i, j, itert, n, k ;
int m = 0;
double v[7];
double *q0, *q, *q1, *qh, *z, *z0, *z1, *zh, *cko, *areao;
double lambda ;
double delta_x, delta_t;
double b , length, qq0, bb, tt, h0, t, g = 9.81, rn=0.023;
FILE *fp;
fp = fopen("data.txt", "r");
if (fp == NULL)
return 0;
while (fscanf(fp, "%lf", &v[m]) != EOF)
m++;
fclose(fp);
b = v[0];
h0 = v[1];
qq0 = v[2];
length = v[3];
delta_x = v[4];
delta_t = v[5];
t = v[6];
bb = length / delta_x;
tt = t / delta_t;
itert = floor(bb);
k = round(tt);
n = itert;
printf("n=%d,k=%d\n", n, k);
printf("b=%f,h0=%f,qq0=%f\n", b, h0, qq0);
printf("length=%f,delta_x=%f,delta_t=%f,t=%f\n", length, delta_x, delta_t, t);
q0 = (double *)malloc((n + 1) * sizeof(double));
q = (double *)malloc((n + 1) * sizeof(double));
q1 = (double *)malloc((n + 1) * sizeof(double));
qh = (double *)malloc((n + 1) * sizeof(double));
z0 = (double *)malloc((n + 1) * sizeof(double));
z = (double *)malloc((n + 1) * sizeof(double));
z1=(double *)malloc((n + 1) * sizeof(double));
zh = (double *)malloc((n + 1) * sizeof(double));
areao = (double *)malloc((n + 1) * sizeof(double));
cko = (double *)malloc((n + 1) * sizeof(double));
for (i = 1; i <= n; i++)
{
z[i] = 1;
q[i] = 1;
}
for (i = 1; i <= n; i++)
{
cko[i] = ck(b, z[i], rn);
areao[i] = area(b, z[i]);
}
q[0] = qq0;
q[n+1] = q[1];
z[0] = h0;
z[n + 1] = z[1];
areao[0] = b * z[0];
areao[n + 1]= areao[1];
for (i = 1; i <= n+1; i++)
{
q0[i] = q[i];
z0[i] = z[i];
}
lambda = delta_t / delta_x;
//First time step is FTBS
for (i = 1; i <= n; i++)
{
zh[i] = z[i] - lambda * (1 / b)*(q[i] - q[i - 1]);
qh[i] = (1 / (g*delta_t*((areao[i] * fabs(q[i])) / (cko[i] * cko[i])) + 1))*(q[i] * (1 - (g*delta_t*((areao[i] * fabs(q[i])) / (cko[i] * cko[i])))) - 2 * lambda*(q[i] / areao[i])*(q[i] - q[i - 1]) + lambda * ((q[i] * q[i]*b - g * areao[i] * areao[i] * areao[i]) / (areao[i] * areao[i])*(z[i] - z[i - 1])) + lambda * ((q[i] * q[i]) / (areao[i] * areao[i]))*z[i] * (areao[i] - areao[i - 1]));
}
for (i = 1; i <= n; i++)
{
printf("%f,%f\n", zh[i], qh[i]);
}
for (i = 1; i <= n; i++)
{
z1[i] = zh[i];
q1[i] = qh[i];
}
for (i = 1; i <= n; i++)
{
areao[i] = area(b, z1[i]);
}
q1[0] = q1[n];
q1[n + 1] = q1[n];
z1[0] = z1[n];
z1[n + 1] = z1[n];
areao[0]= b * z1[0];
areao[n+1] = areao[n];
for (j = 0; j < k; j++)
{
for (i = 1; i <= n; i++)
{
zh[i] = z[i] - lambda * (1 / b)*(q1[i + 1] - q1[i - 1]);
qh[i] = (1/(g*delta_t*((areao[i] * fabs(q1[i])) / (cko[i] * cko[i])) + 1))*(q[i] * (1 - (g*delta_t*((areao[i] * fabs(q1[i])) / (cko[i] * cko[i])))) - 2 * lambda*(q1[i] / areao[i])*(q1[i+1] - q1[i - 1]) + lambda * ((q1[i] * q1[i] * b - g * areao[i] * areao[i] * areao[i]) / (areao[i] * areao[i])*(z1[i+1]-z1[i-1])) + lambda * ((q1[i] * q1[i]) / (areao[i] * areao[i]))*z1[i] * (areao[i+1] - areao[i - 1]));
}
for (i = 0; i <= n + 1; i++)
{
z[i] = z1[i];
q[i] = q1[i];
}
for (i = 1; i <= n; i++)
{
z1[i] = zh[i];
q1[i] = qh[i];
}
q1[0] = q1[n];
q1[n + 1] = q1[n];
z1[0] = z1[n];
z1[n + 1] = z1[n];
for (i = 1; i <= n; i++)
{
printf("%d %f,%f,%f,%f\n", i, z0[i], z1[i], q0[i], q1[i]);
}
}
}
结果