nonreflection boundary conditions

无反射边界条件数值计算

一．问题描述

一维欧拉方程：

0U F

t x

∂∂+=∂∂ (1)

()2,u U u F u p e e p u ρρρρ⎛⎫

⎛⎫

⎪ ⎪=+ ⎪ ⎪ ⎪ ⎪+⎝⎭⎝⎭

(2)

通过差分方法在计算域a x b ≤≤上求解方程(1)时，需要知道a x b ≤≤上的初始条件和,x a x b ==的边界条件。无反射边界条件即讨论如何处理边界条件的问题。 通过线性变换，方程（1）可以写作： 0U U

A t x

∂∂+=∂∂ (3)

其中：

20,00

u

c p U u A u s s u ρ

ρρρ⎛⎫

⎛⎫

⎪ ⎪

⎪

== ⎪ ⎪ ⎪ ⎪⎝⎭ ⎪⎝

⎭

(4)

求矩阵A 的特征值i λ，及其相应的左右特征向量,i i I r 。

i i i i i i

I A I Ar r λλ== (5)

自然可以得到对角变换： 1SAS -=Λ

(10)

其中：

()11

121232

33,,,,I S I S r r r I λλλ-⎛⎫⎛⎫

⎪ ⎪==Λ= ⎪ ⎪ ⎪ ⎪⎝⎭⎝

⎭

(11)

(10)带入(3)得到： 0U U

S

S t x ∂∂+Λ=∂∂ (12) 0i i i U U I I t x λ∂∂+=∂∂ (13)

令：

i i dV I dU =

(14)

所以：

0i i i V V

t x

λ∂∂+=∂∂ (15)

式(15)描述了一系列特征速度为i λ的波动方程。对于方程(3)变换以后得到(15),其中：

()1231232123, , ,,, 0,0,1, ,,, , u c u u c

p p I c I I c sc sc V p cu V p c V p cu

λλλρρρρρ=-==+⎛⎫⎛

⎫=--== ⎪ ⎪⎝⎭⎝

⎭=-=-=+(16)

现在来看如何根据(15)来给定计算过程中的边界条件，以计算域的左侧边界x a =为例。当波速0i λ>，表示波由计算域外传入，而实际上若x a =点处并没有实际的源存在，则应当有

0i V x ∂=∂，即0i V

t

∂=∂，这就保证不会在无源有计算所引进的误差而使有i V 传入计算域，若x a =点本身就有源，则直接根据给定的源确定边界条件即可；当0i λ<，表示波右计算域内往外传播，这个时候就应当由前差的方法，由x a =右侧点n t t =时刻的信息来得到x a =点在n t t t =+∆时刻的边界条件。一下给出具体的取值格式： 离散方程组： 10i i i i i dp du

c dt dt ρη-+= (17) 220i i

i i dp d c dt dt ρη-+= (18)

30i i

i i i dp du c dt dt

ρη++= (19)

其中j η分别为相应j V 的空间导数，当在边界处无源时，应当安一下方法取值。若波由外传入，则取0j η=，若波由计算域传出，则有：

,1,j ,,1j []/ 0[]/ 0

ji j j i j i ji ji j i j i V V dx V V dx ηλληλλ+-=-<=-> (20)

特征值方法处理无反射边界条件可以很好的帮助理解，但在做高维问题是会比较

困难。

二.作业：

一维直管道（如下图所示），计算域为 1010≤≤-x ，在下游出口处有一声源为：

)]1(

sin[111'''t Ma x

p u +-⎪

⎪⎪⎭

⎫ ⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛ωερ， 510-=ε， πω6.0= 主控方程为一维线化Euler 方程，计算t=1, 3, 5, 10, 20, 40，60，80等时刻管内声场分布，并与解析解对比。

三．数值程序

#include #include #include #include #include #include

using namespace std;

const int N=1001;

const double gamma=1.4; const double EPS=0.00001; const double Ma=0.4; const double XR=10.0; const double XL=-10.0;

const double PI=3.1415926538; int main() {

int i,step;

double dx,dt,t,w,c0,rho0,p0,u0;

double x[N],rho[N],p[N],u[N],lambda[3],VV[N][3]; w=0.6*PI; dt=0.0001;

dx=(XR-XL)/(N-1);

for(i=0;i<=N-1;i++)

x[i]=XL+i*dx;

c0=1.0;

rho0=1.0;

u0=Ma;

p0=1.0/gamma;

for(i=0;i<=N-1;i++)

{

rho[i]=0.0;

u[i]=Ma;

p[i]=0.0;

VV[i][0]=0.0;

VV[i][1]=0.0;

VV[i][2]=0.0;

}

step=0;

lambda[0]=Ma-1;

lambda[1]=Ma;

lambda[2]=Ma+1;

t=0.0;

do

{

step++;

//------right boundary----

rho[N-1]=EPS*sin(w*(XR/(1.0-Ma)+t));

u[N-1]=-EPS*sin(w*(XR/(1.0-Ma)+t));

p[N-1]=EPS*sin(w*(XR/(1.0-Ma)+t));

VV[N-1][0]=p[N-1]-rho0*c0*u[N-1];

VV[N-1][1]=p[N-1]-c0*c0*rho[N-1];

VV[N-1][2]=p[N-1]+rho0*c0*u[N-1];

//------left boundary----

VV[0][0]=VV[0][0]-lambda[0]*dt*(VV[1][0]-VV[0][0])/dx;

VV[0][1]=VV[1][1];

VV[0][2]=VV[1][2];

//-----inner flow-------

for(i=N-2;i>=1;i--)

{

VV[i][0]=VV[i][0]-0.5*lambda[0]*dt*(VV[i+1][0]-VV[i-1][0])/dx;

VV[i][1]=VV[i][1]-0.5*lambda[1]*dt*(VV[i+1][1]-VV[i-1][1])/dx;