1、谱元方法简介,主讲：秦国良 电话：82663537 邮箱：,2019年6月6日,Tel：82663537 e-mail：,2,1、正交函数系与谱近似,1.1 正交函数系与正交多项式 1.2 函数的Fourier展开 1.3 Chebyshev谱逼近(离散函数的Fourier展开) 1.4 插值函数的导数 1.5 二维函数的Chebyshev谱近似,2019年6月6日,Tel：82663537 e-mail：,3,1.1 正交函数系与正交多项式,设函数系j(x)，如满足条件，则称函数系j(x)，在区间a,b上关于权函数w(x) 正交， 当j(x)是代数多项式时，称为正交多项式。如Chebyshe
2、v多项式，,2019年6月6日,Tel：82663537 e-mail：,4,Chebyshev多项式的正交性,2019年6月6日,Tel：82663537 e-mail：,5,1.2 函数的Fourier展开,函数f(x)在区间a,b上满足Dirichlet条件，且加权平方可积，即对于权函数w(x)，有，则f(x)可在a,b上以w(x)为权函数，按正交多项式n(x)展成广义Fourier级数，设j(x)为在节点xk,k=0,1,N上的正交函数系，权为wk0,k=0,1,N，设f(x)为在节点系xk, k = 0,1,N上取值的已知函数，则函数f(x)的Fourier展开定义为，,2019年6
3、月6日,Tel：82663537 e-mail：,6,1.3 Chebyshev谱逼近(离散函数的Fourier展开),取节点系xk, k = 0,1,N为N阶Chebyshev多项式的极值点(称Gauss-Lobatto点)，即，权系数则上面的展开式变为，此式亦可看作f(x)在配置点上的插值, g(x)为插值函数。,2019年6月6日,Tel：82663537 e-mail：,7,1.4 插值函数的导数,对 微分得：插值函数在配置点的导数可表示为，,2019年6月6日,Tel：82663537 e-mail：,8,1.5 二维函数的Chebyshev谱近似,对二维函数u(,)，在标准正方形单
4、元内，定义,方向节点系j，j=0,1,Nx和k ，k=0,1,Ny，,2019年6月6日,Tel：82663537 e-mail：,9,2 谱方法求解微分方程,2.1 解微分方程的加权残量法MWR（Method of Weight Residuals) 2.2 谱方法求解微分方程 2.3 小结,2019年6月6日,Tel：82663537 e-mail：,10,2.1 解微分方程的加权残量法MWR（Method of Weight Residuals),基本思想 MWR的步骤 MWR的基本方法,2019年6月6日,Tel：82663537 e-mail：,11,加权残量法基本思想(1),设微分
5、方程边值问题为，加权余量法，是求微分方程形式如下的近似解。a=a0,a1,aNT为待定向量，0, 1, N为为基函数，且是一组线性无关的函数,2019年6月6日,Tel：82663537 e-mail：,12,加权残量法基本思想(2),“残量”定义如下，显然只有当试函数为边值问题的精确解时，余量(2.31)才为零。 加权残量法 适当选取待定向量a=a0,a1,aNT，使得“残量”极小。通常某一权函数系wi(x),i=0,1,N，使得权函数与“残量”正交，来确定待定向量。,2019年6月6日,Tel：82663537 e-mail：,13,MWR的步骤,选取一组满足要求的基函数构造试函数 并满足
6、所有边界条件；选取一组权函数运用MWR准则， 得到 的代数方程组； 解上述代数方程组，确定待定参数。,2019年6月6日,Tel：82663537 e-mail：,14,MWR的基本方法,按权函数的不同有以下几种基本方法 Galerkin法(Galerkin Method)配置法(Collocation Method)Tan方法(Tan Method)类似于Galerkin法，Tan方法中，试函数不需要满足边界条件，增加边界约束系数来满足边界条件,2019年6月6日,Tel：82663537 e-mail：,15,2.2 谱方法求解微分方程,以非线性Burger方程为例来说明谱方法的离散过程,
7、并区别Galerkin法,Tau方法和Collocation方法。Burger方程如下:,2019年6月6日,Tel：82663537 e-mail：,16,Chebyshev Galerkin方法 (Galerkin Method,方程的近似解表达为， 由MWR准则得,2019年6月6日,Tel：82663537 e-mail：,17,Chebyshev Tau方法(Tan Method,方程的近似解表达为， 由MWR准则得, 边界约束,2019年6月6日,Tel：82663537 e-mail：,18,Chebyshev Collocation方法(Collocation Method),
8、选取节点系为N阶Chebyshev多项式的极值点, 方程的近似解用下式表示，由Collocation方程得,2019年6月6日,Tel：82663537 e-mail：,19,2.3 小结,一种求解偏微分方程的高阶精度数值方法。 属于求解偏微分方程加权残量法的特例。 谱方法使用高阶正交多项式作为展开函数 。 谱方法最受人青睐的优越性在于它具有“无穷阶”的收敛速度，其确切含义为，若原微分方程的解无限可微，则由适当的谱方法所得到的近似解对原问题的收敛速度比1/N的任何幂都更快，这里N是所取基函数的个数。 Kreiss和Oliger, Orszag Gottlieb和Orszag,2019年6月6日
9、,Tel：82663537 e-mail：,20,3 谱元方法求解微分方程,回顾谱方法 使用加权残量法 谱元方法基本思想 常微分方程边值问题的Galerkin变分原理 偏微分方程的Galerkin变分原理 Galerkin逼近解 常微分方程元素矩阵的形成 偏微分方程元素矩阵的形成 极坐标系下环形区域的谱元方法 总刚度矩阵的形成,2019年6月6日,Tel：82663537 e-mail：,21,3.1 回顾谱方法 使用加权残量法,选取一组满足要求的基函数,构造试函数，满足所有的边界条件,选取一组权函数,利用MWR准则，得到ai 的代数方程,求解上述代数方程组，确定待定参数,2019年6月6日,
10、Tel：82663537 e-mail：,22,3.2 谱元方法基本思想,使用有限元方法的思想 ，将求解域分成若干子域 采用谱方法中谱近似技术代替有限元中的插值函数 采用有限元中等参单元,亦可逼近复杂求解边界,2019年6月6日,Tel：82663537 e-mail：,23,3.3 常微分方程边值问题的Galerkin变分原理,考察二点边值问题上述问题的Galerkin变分问题是：求 ，使得，其中,2019年6月6日,Tel：82663537 e-mail：,24,3.4 偏微分方程的Galerkin变分原理,考察二维椭圆边值问题上述问题的Galerkin变分问题是：求 ，使得，,2019年
11、6月6日,Tel：82663537 e-mail：,25,3.5 Galerkin逼近解,设Vh是V的有限维子空间，当h-0时，Vh的维数无限增加，直到充满V为止，那么Galerkin逼近解 ，使得，设 是Vh的基函数系，那么由于 ，故有，其中得由于 的任意性，得，即得线性代数方程组，可以证明此代数方程组的解唯一存在，即可解出Galerkin逼近解。,2019年6月6日,Tel：82663537 e-mail：,26,3.6 常微分方程元素矩阵的形成,考察以下微分方程用谱方法求解此微分方程,先将求解区域分成若干单元,设单元总数为N,考察第i个单元ei,将ei转化为标准求解单元,设x为总体坐标,
12、为标准求解单元内的局部坐标,并设,2019年6月6日,Tel：82663537 e-mail：,27,插值点和插值函数,在第i个元素内选取Nx+1个点作为u(x)和f(x)的插值点,本文选取阶Chebyshev多项式的极值点作为插值点,即 ,u(x)和f(x)可以表达为，,2019年6月6日,Tel：82663537 e-mail：,28,元素矩阵形成,对一维边值问题,在ei单元内线性方程式变为 其中，上面线性方程组变为，,2019年6月6日,Tel：82663537 e-mail：,29,的推导过程,2019年6月6日,Tel：82663537 e-mail：,30,3.7 偏微分方程元素矩
13、阵的形成,具体考察Helmholtz方程，其中， 作以下坐标变换在标准正方形单元内，使用伪谱Chebyshev谱逼近u(,)，则元素内线性方程组式变为，,2019年6月6日,Tel：82663537 e-mail：,31,总体坐标系与局部坐标系,2019年6月6日,Tel：82663537 e-mail：,32,元素矩阵,2019年6月6日,Tel：82663537 e-mail：,33,3.8 极坐标系下环形区域的谱元方法,在极坐标系下，Helmholtz方程表达为，将极坐标系下环形段单元，变换为直角坐标系下的正方形单元，坐标变换关系如下，元素矩阵的形式为，,2019年6月6日,Tel：82
14、663537 e-mail：,34,极坐标坐标系与局部坐标系,2019年6月6日,Tel：82663537 e-mail：,35,元素矩阵,式中各变量的表达式如下，,2019年6月6日,Tel：82663537 e-mail：,36,3.9 总刚度矩阵的形成,节点影响元素集，节点所在的元素的并集 影响节点，对某点，其影响元素内的所有节点 相互影响元素集，两个节点和影响元素之交集 对节点和，总体矩阵元素为，由于 在一个元素e内计算，故称它为单元矩阵，上式表明矩阵有贡献的只是()影响元素集所有单元矩阵相应的元素。,2019年6月6日,Tel：82663537 e-mail：,37,4 复杂区域的谱
15、元方法,采用谱逼近的等参数单元 总体坐标系下导数的计算 坐标变换矩阵及变换行列式 局部坐标系下偏导数的计算 总体坐标系下偏导数的计算 复杂求解区域的谱元方法,2019年6月6日,Tel：82663537 e-mail：,38,总体坐标与局部坐标系,2019年6月6日,Tel：82663537 e-mail：,39,4.1采用谱逼近的等参数单元,总体坐标与局部坐标(如图所示) 等参单元，插值函数表达式和坐标变换表达式一致 函数的插值表达为，坐标变换的表达式为，导数表达式,2019年6月6日,Tel：82663537 e-mail：,40,4.2 总体坐标系下导数的计算,坐标变换矩阵(Jacobi
16、矩阵)及变换行列式(Jacobi行列式),微元面积、微元弧长(为常数)时， (为常数)时,2019年6月6日,Tel：82663537 e-mail：,41,局部坐标系下偏导数的计算,由插值式对求偏导得，上式在(m,n)处的值为，比照函数得，则局部坐标系下的偏导数表达为，,2019年6月6日,Tel：82663537 e-mail：,42,总体坐标系下偏导数的计算,由插值函数得，由Jacobi矩阵得总体坐标系下，偏导数计算式为，,2019年6月6日,Tel：82663537 e-mail：,43,4.3 复杂求解区域的谱元方法,考察Helmholtz方程则元素矩阵为，,2019年6月6日,Te
17、l：82663537 e-mail：,44,元素矩阵各变量表达式,2019年6月6日,Tel：82663537 e-mail：,45,5 非定常不可压缩流动与自然对流换热问题求解,不可压缩流动控制方程 封闭空腔内自然对流换热基本控制方程 解非定常不可压缩Navier-Stokes方程的时间分裂法 时间分裂法求解自然对流问题 时间分裂法求解Navier-Stokes方程的误差分析,2019年6月6日,Tel：82663537 e-mail：,46,5.1 非定常不可压缩流动控制方程,非定常不可压缩流动无因次Navier-Stokes方程为，以顶盖移动的方腔驱动流为例，定义无因次参量如下，,201
18、9年6月6日,Tel：82663537 e-mail：,47,5.2 封闭空腔内自然对流换热基本控制方程,封闭空腔内自然对流换热基本控制方程,2019年6月6日,Tel：82663537 e-mail：,48,无因次量定义,对正方形封闭空腔，无因次参数的定义如下，对空气来说，,2019年6月6日,Tel：82663537 e-mail：,49,5.3 解非定常不可压缩Navier-Stokes方程的时间分裂法,Navier-Stokes可写成如下形式，对时间积分得，设 为三步中间值，则分三步求解的过程为：非线性步、压力步、粘性步,2019年6月6日,Tel：82663537 e-mail：,5
19、0,非线性步,非线性步使用显式4阶Runge-Kutta方法，替代传统时间分裂法的Adams-Bashforth方法，使本时间步仅与上一时间步有关，此步中代入初始条件，不使用边界条件。,2019年6月6日,Tel：82663537 e-mail：,51,压力步,假设 满足不可压缩条件，且在边界上法向上的投影为零，即，得,此步为解带有第二类边界条件的压力Poisson方程，解出压力Poisson方程后，再由上式求解 。,2019年6月6日,Tel：82663537 e-mail：,52,粘性步,利用Crank-Nicolson格式离散粘性步，得Helmholtz方程，施以固体壁面边界条件即可求解
20、,2019年6月6日,Tel：82663537 e-mail：,53,5.4 时间分裂法求解自然对流问题,求解自然对流问题的过程同Navier-Stokes方程，也是将一个时间步分为三步求解，所不同的是在一个时间步内除了解Navier-Stokes外，还需求解能量方程，能量方程中没有压力项，因此在第二步压力Poisson方程求解中与能量方程无关。,2019年6月6日,Tel：82663537 e-mail：,54,非线性步,非线性步仍使用显式4阶Runge-Kutta方法，,2019年6月6日,Tel：82663537 e-mail：,55,压力步,为解带有第二类边界条件的压力Poisson方
21、程，解出压力Poisson方程后，再由(5.34)式求解 。,2019年6月6日,Tel：82663537 e-mail：,56,粘性步,利用Crank-Nicolson格式离散粘性步，得Helmholtz方程，施以固体壁面边界条件即可求解。,2019年6月6日,Tel：82663537 e-mail：,57,5.5 时间分裂法求解Navier-Stokes方程的误差分析,利用时间分裂法求解Navier-Stokes方程基于以下两个假设，满足不可压缩约束，即满足连续性方程，在边界的法线上的投影为零，高精度压力边界条件描述，压力Poisson方程变为，,2019年6月6日,Tel：8266353
22、7 e-mail：,58,6 面向对象谱元方法程序设计,面向对象程序设计思想 面向对象谱元方法程序设计 节点类 标准单元类 单元类 整体区域类 谱元方法程序设计中的几个问题 插值函数积分的存储 节点的编码及总体矩阵的形成 Navier-Stokes方程和能量方程的求解过程,2019年6月6日,Tel：82663537 e-mail：,59,6.1 面向对象程序设计思想,按人们通常的思维方式建立模型，设计尽可能自然地表现求解方法的软件。为此，必须建立直接表现问题中事物与事物相互联系的概念，建立适合人们思维方式的描述方法；在面向对象的方法中主要有以下概念：对象、消息传递、类和继承、方法等。 对象和
23、消息传递是表示事物与事物间的关系； 类和继承是适应人们思维方式对事物的描述方法； 方法是作用在对象上的各种操作。,2019年6月6日,Tel：82663537 e-mail：,60,6.2 面向对象谱元方法程序设计,面向对象的程序设计包括面向对象分析、面向对象设计和面向对象实现对谱元方法来说，从谱元方法分析问题的思路出发，找出描述谱元方法的对象类（数据和方法的结合体），建立对象类的关系，并集成这些对象。设计节点类、标准单元、单元类、总体区域类。这些类中用属性来描述类的特性，用方法来实现某些操作。有些类中的属性可能又需要用对象类来描述，方法的作用对象也可能是某一类的实体。,2019年6月6日,T
24、el：82663537 e-mail：,61,节点类,主要的数据及方法包含节点的属性，如节点的编号、节点的位置坐标、节点的边界属性、微分方程右端项及边界条件值、所要求解的变量及其导数值。 节点类定义， class CGridPoint protected: /节点的属性标准节点，编码，坐标，边界属性，方程右端项；压力，速度及其导数； public: /节点类方法属性设置，边界设置，压力和速度设置；属性取值，边界属性取值； ,2019年6月6日,Tel：82663537 e-mail：,62,标准单元类,包括了在各求解单元中的常量。类中的数据及方法主要有：单元内网格的划分(含两方向的节点数及总节
25、点数)、标准插值节点、插值函数的导数矩阵及有关标准插值函数的积分。 标准单元类定义， class CStdElement protected: /标准单元类的属性网格属性，标准插值节点，插值函数的导数，插值函数积分； public: /标准单元类方法标准单元类的构造，Chebyshev多项式及其导数；插值函数的积分，属性设置及取值； ,2019年6月6日,Tel：82663537 e-mail：,63,单元类,主要数据及方法包含元素序号、单元内节点划分(和标准单元内划分相对应)、标准单元类、节点类和总体坐标与局部坐标的变换关系。 单元类定义， class CElement protected:
26、 /单元类的属性单元及网格属性，标准单元类，坐标变换，单元矩阵； public: /单元类方法单元构造，标准单元类的构造，节点类构造；元素矩阵的形成，属性设置及取值； ,2019年6月6日,Tel：82663537 e-mail：,64,总体区域类,主要数据及方法包含区域内单元数、总节点数、节点类、标准单元类、总体矩阵等。 总体区域类定义， class CDomain protected: /整体区域类的属性区域的属性，单元及网格属性，单元类，标准单元类；总体矩阵，对角元位置，边界节点序号； public: /整体区域类方法单元构造，标准单元类的构造，节点类构造；元素矩阵的形成，属性设置及取值
27、； ,2019年6月6日,Tel：82663537 e-mail：,65,6.2 谱元方法程序设计中的几个问题,插值函数积分的存储 节点的编码及总体矩阵的形成 对称带状矩阵的一维存储 元素的编码 影响元素和影响节点 对角元的位置及一维数组的长度 总刚度矩阵 第一类边界条件的处理,2019年6月6日,Tel：82663537 e-mail：,66,插值函数积分的存储,在形成元素矩阵时，要计算如下积分，可以看出，在(i,j,k)变化时，不必占用(Nx+1)3个存储单元，这样既节约了存储空间又节省了计算时间。Bijk所占的空间可由下式计算，Bijk的检索按如下方式进行，先将i,j,k按由大到小排列，
28、设ijk，则Bijk在数组中的位置为，,2019年6月6日,Tel：82663537 e-mail：,67,6.3 Navier-Stokes方程和能量方程的求解过程,构造标准单元stdElement，求出对各单元不变的量。 构造求解域pGeometry，对求解域进行网格划分，同时构造节点类pDotNG，构造各求解单元pElementNe，对节点类中U、V、赋初值，赋零值，在热边界上赋值1，在冷边界上赋0，其余两边为第二类边界条件，在形成单元矩阵时自动满足。 利用Runge-Kutta法求解 、 、 ，将求解结果存于 、 、 解压力Poisson方程，第二类边界条件，解出压力后更新 、 、 求
29、解关于 、 、 的Helmholtz方程，边界条件第一或第二类。解出Helmholtz方程后更新 、 、 ，即完成了一个时间步的求解。,2019年6月6日,Tel：82663537 e-mail：,68,7 计算实例,常微分方程求解 Helmholtz和Poisson的求解 简单边界Helmholtz方程求解 复杂边界Poisson的求解 移动顶盖方腔驱动非定常流动 封闭方腔中的自然对流,2019年6月6日,Tel：82663537 e-mail：,69,7.1 常微分方程求解,【例1】用谱方法求解下面方程该问题的精确解的表达式为： 将求解域0,-1等分为10个求解单元,单元内使用5阶Cheb
30、yshev多项式进行谱逼近。数值计算结果与精确解的比较如下,比较表明谱元方法具有相当高的精度。,2019年6月6日,Tel：82663537 e-mail：,70,例1 数值计算结果与精确解的比较如下,2019年6月6日,Tel：82663537 e-mail：,71,【例2】,【例2】用谱元方法求解下面方程该问题的精确解的表达式为： 将求解域0,-1等分为4个求解单元,单元内使用5阶Chebyshev多项式进行谱逼近。数值计算结果与精确解的比较如下,比较表明谱元方法具有相当高的精度,2019年6月6日,Tel：82663537 e-mail：,72,例2 数值计算结果与精确解的比较如下,20
31、19年6月6日,Tel：82663537 e-mail：,73,7.2 简单边界Helmholtz方程求解,具体考察Helmholtz方程，精确解为，采用矩形单元采用曲边四边形单元,2019年6月6日,Tel：82663537 e-mail：,74,两种单元求解结果比较(简单边界),表7-1(采用矩形单元的求解结果)表7-2(采用曲边四边形单元的求解结果),2019年6月6日,Tel：82663537 e-mail：,75,复杂边界Poisson的求解,考察方程边界条件由精确解给出采用环形段单元 采用曲边四边形单元,2019年6月6日,Tel：82663537 e-mail：,76,两种单元求
32、解结果比较(复杂边界),表7-3(采用环形段单元的求解结果)表7-4(采用曲边四边形单元的求解结果),2019年6月6日,Tel：82663537 e-mail：,77,7.3 移动顶盖方腔驱动非定常流动,是计算流体力学典型数值实验模型，在很多情况下都是用它来检验数值方法的有效性；除移动顶盖外，初始速度场为零；移动顶盖上速度u=1,v=0，其余边界速度为零；取9个求解单元，使用8阶Chebyshev多项式插值；计算了Re=100的情况，取两种时间步长T1=0.004，T2=0.005，计算2000和1600步，各需12和9.6小时,2019年6月6日,Tel：82663537 e-mail：,
33、78,移动顶盖方腔驱动流物理模型,2019年6月6日,Tel：82663537 e-mail：,79,竖直中心线上的速度U, 水平中心线上的速度V,图7-3 竖直中心线上的速度U,图7-4 水平中心线上的速度V,2019年6月6日,Tel：82663537 e-mail：,80,方腔内中心点的速度U 和 速度V随时间的变化,图7-5 方腔内中心点的速度U,图7-6 方腔内中心点的速度V,2019年6月6日,Tel：82663537 e-mail：,81,封闭方腔中的自然对流,计算中网格单元数为3x3，网格内部的插值函数采用8阶Chebyshev多项式；计算的空气在正方形封闭空腔内的自然对流，空
34、气的Pr数取为0.71Ra=103, 时间步长取T=3x10-4，计算1500时间步；耗时12小时 Ra=104 , 时间步长取T=3x10-4，计算1500时间步；耗时12小时 Ra=105 , 时间步长取T=4x10-5，计算12500时间步；耗时90小时,2019年6月6日,Tel：82663537 e-mail：,82,封闭方腔中的自然对流物理模型,2019年6月6日,Tel：82663537 e-mail：,83,数值模拟结果与基准解的比较,数值模拟结果与基准解的比较,2019年6月6日,Tel：82663537 e-mail：,84,竖直中心线u速度分布、水平中心线上v速度分布,图
35、7-13 竖直中心线u速度分布,图7-14水平中心线上v速度分布,2019年6月6日,Tel：82663537 e-mail：,85,水平中心线上的温度分布,图7-15水平中心线上的温度分布,2019年6月6日,Tel：82663537 e-mail：,86,Ra=103速度矢量、温度场,2019年6月6日,Tel：82663537 e-mail：,87,Ra=104的速度矢量和温度场,2019年6月6日,Tel：82663537 e-mail：,88,Ra=105的速度矢量和温度场,总结,函数逼近 谱元方法和有限元方法 时间分裂法解N-S方程 发展谱元方法的几个问题,2019年6月6日,Te
36、l：82663537 e-mail：,89,2019年6月6日,Tel：82663537 e-mail：,90,了解谱元方法,谱元,谱方法,有限元方法,Spectral,Element,2019年6月6日,Tel：82663537 e-mail：,91,有限元类方法求解过程,方程变分,区域划分,定基函数,单元分析,合成总刚,边界条件,问题的解,原始方程,2019年6月6日,Tel：82663537 e-mail：,92,谱元方法的基本原理,2019年6月6日,Tel：82663537 e-mail：,93,时间分裂法的应用,2019年6月6日,Tel：82663537 e-mail：,94,谱元方法的若干问题,时间方向的精度问题（高Re数问题）,复杂区域的进一步适应性,计算速度问题,圆形区域问题,一般区域问题,极坐标或者柱坐标系,等参谱元方法,提高时间分 裂法的阶数,提高压力边 界条件阶数,
