谱元方法课程

谱元方法课程第1页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日21、正交函数系与谱近似1.1正交函数系与正交多项式1.2函数的Fourier展开1.3Chebyshev谱逼近(离散函数的Fourier展开)1.4插值函数的导数1.5二维函数的Chebyshev谱近似第2页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日31.1正交函数系与正交多项式设函数系{j(x)}，如满足条件，则称函数系{j(x)}，在区间[a,b]上关于权函数w(x)正交，当j(x)是代数多项式时，称为正交多项式。如Chebyshev多项式，

第3页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日4Chebyshev多项式的正交性第4页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日51.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}上的正交函数系，权为{wk>0,k=0,1,…,N}，设f(x)为在节点系{xk,k=0,1,…,N}上取值的已知函数，则函数f(x)的Fourier展开定义为，第5页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日61.3Chebyshev谱逼近(离散函数的Fourier展开)取节点系{xk,k=0,1,…,N}为N阶Chebyshev多项式的极值点(称Gauss-Lobatto点)，即，权系数则上面的展开式变为，此式亦可看作f(x)在配置点上的插值,g(x)为插值函数。第6页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日71.4插值函数的导数对微分得：插值函数在配置点的导数可表示为，第7页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日81.5二维函数的Chebyshev谱近似对二维函数u(ξ,η)，在标准正方形单元内，定义ξ,η方向节点系{ξj，j=0,1,…,Nx}和{ηk，k=0,1,…,Ny}，第8页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日92谱方法求解微分方程2.1解微分方程的加权残量法MWR

（MethodofWeightResiduals)2.2谱方法求解微分方程

2.3小结第9页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日102.1解微分方程的加权残量法MWR（MethodofWeightResiduals)基本思想MWR的步骤MWR的基本方法第10页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日11加权残量法基本思想(1)设微分方程边值问题为，加权余量法，是求微分方程形式如下的近似解。a=[a0,a1,…,aN]T为待定向量，0,1,…,N为为基函数，且是一组线性无关的函数第11页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日12加权残量法基本思想(2)“残量”定义如下，显然只有当试函数为边值问题的精确解时，余量(2.31)才为零。加权残量法适当选取待定向量a=[a0,a1,…,aN]T，使得“残量”极小。通常某一权函数系{wi(x),i=0,1,…,N}，使得权函数与“残量”正交，来确定待定向量。

第12页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日13MWR的步骤①选取一组满足要求的基函数②构造试函数并满足所有边界条件；③选取一组权函数④运用MWR准则，得到的代数方程组；⑤解上述代数方程组，确定待定参数。第13页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日14MWR的基本方法按权函数的不同有以下几种基本方法Galerkin法(GalerkinMethod)配置法(CollocationMethod)Tan方法(TanMethod)

类似于Galerkin法，Tan方法中，试函数不需要满足边界条件，增加边界约束系数来满足边界条件

第14页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日152.2谱方法求解微分方程

以非线性Burger方程为例来说明谱方法的离散过程,并区别Galerkin法,Tau方法和Collocation方法。Burger方程如下:

第15页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日16ChebyshevGalerkin方法

(GalerkinMethod方程的近似解表达为，

由MWR准则得,

第16页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日17ChebyshevTau方法(TanMethod方程的近似解表达为，

由MWR准则得,

边界约束,第17页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日18ChebyshevCollocation方法(CollocationMethod)选取节点系为N阶Chebyshev多项式的极值点,

方程的近似解用下式表示，由Collocation方程得,第18页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日192.3小结一种求解偏微分方程的高阶精度数值方法。属于求解偏微分方程加权残量法的特例。谱方法使用高阶正交多项式作为展开函数。谱方法最受人青睐的优越性在于它具有“无穷阶”的收敛速度，其确切含义为，若原微分方程的解无限可微，则由适当的谱方法所得到的近似解对原问题的收敛速度比1/N的任何幂都更快，这里N是所取基函数的个数。Kreiss和Oliger,OrszagGottlieb和Orszag第19页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日203谱元方法求解微分方程回顾谱方法———使用加权残量法谱元方法基本思想常微分方程边值问题的Galerkin变分原理偏微分方程的Galerkin变分原理Galerkin逼近解常微分方程元素矩阵的形成偏微分方程元素矩阵的形成极坐标系下环形区域的谱元方法总刚度矩阵的形成第20页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日213.1回顾谱方法———使用加权残量法

选取一组满足要求的基函数

构造试函数，满足所有的边界条件

选取一组权函数

利用MWR准则，得到ai的代数方程

求解上述代数方程组，确定待定参数第21页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日223.2谱元方法基本思想使用有限元方法的思想，将求解域分成若干子域

采用谱方法中谱近似技术代替有限元中的插值函数

采用有限元中等参单元,亦可逼近复杂求解边界

基函数的选用在一个单元上选取在整个区域上选取计算精度不是太高比较高单元的选择三角形，四边形等四边形提高精度方法网格加密不改变插值次数增加基函数多项式的次数第22页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日233.3常微分方程边值问题的Galerkin变分原理考察二点边值问题上述问题的Galerkin变分问题是：求，使得，其中第23页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日243.4偏微分方程的Galerkin变分原理考察二维椭圆边值问题上述问题的Galerkin变分问题是：求，使得，第24页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日253.5Galerkin逼近解设Vh是V的有限维子空间，当h->0时，Vh的维数无限增加，直到充满V为止，那么Galerkin逼近解，使得，设是Vh的基函数系，那么由于，故有，其中得由于的任意性，得，即得线性代数方程组，可以证明此代数方程组的解唯一存在，即可解出Galerkin逼近解。第25页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日263.6常微分方程元素矩阵的形成

考察以下微分方程用谱方法求解此微分方程,先将求解区域分成若干单元,设单元总数为N,考察第i个单元ei,将ei转化为标准求解单元,设x为总体坐标,为标准求解单元内的局部坐标,并设

第26页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日27插值点和插值函数在第i个元素内选取Nx+1个点作为u(x)和f(x)的插值点,本文选取阶Chebyshev多项式的极值点作为插值点,即

,u(x)和f(x)可以表达为，第27页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日28元素矩阵形成对一维边值问题,在ei单元内线性方程式变为

其中，上面线性方程组变为，第28页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日29的推导过程第29页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日303.7偏微分方程元素矩阵的形成具体考察Helmholtz方程，其中，作以下坐标变换在标准正方形单元内，使用伪谱Chebyshev谱逼近u(ξ,η)，则元素内线性方程组式变为，第30页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日31总体坐标系与局部坐标系第31页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日32元素矩阵第32页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日333.8极坐标系下环形区域的谱元方法在极坐标系下，Helmholtz方程表达为，将极坐标系下环形段单元，变换为直角坐标系下的正方形单元，坐标变换关系如下，元素矩阵的形式为，第33页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日34极坐标坐标系与局部坐标系第34页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日35元素矩阵式中各变量的表达式如下，第35页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日363.9总刚度矩阵的形成节点影响元素集，节点所在的元素的并集影响节点，对某点，其影响元素内的所有节点相互影响元素集，两个节点和影响元素之交集对节点和，总体矩阵元素为，由于在一个元素e内计算，故称它为单元矩阵，上式表明矩阵有贡献的只是()影响元素集所有单元矩阵相应的元素。第36页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日374复杂区域的谱元方法采用谱逼近的等参数单元总体坐标系下导数的计算坐标变换矩阵及变换行列式局部坐标系下偏导数的计算总体坐标系下偏导数的计算复杂求解区域的谱元方法第37页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日38总体坐标与局部坐标系第38页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日394.1采用谱逼近的等参数单元总体坐标与局部坐标(如图所示)等参单元，插值函数表达式和坐标变换表达式一致函数的插值表达为，坐标变换的表达式为，导数表达式第39页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日404.2总体坐标系下导数的计算坐标变换矩阵(Jacobi矩阵)及变换行列式(Jacobi行列式)微元面积、微元弧长

(为常数)时，(为常数)时第40页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日41局部坐标系下偏导数的计算由插值式对ξ求偏导得，上式在(ξm,ηn)处的值为，比照函数得，则局部坐标系下的偏导数表达为，第41页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日42总体坐标系下偏导数的计算由插值函数得，由Jacobi矩阵得总体坐标系下，偏导数计算式为，第42页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日434.3复杂求解区域的谱元方法考察Helmholtz方程则元素矩阵为，第43页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日44元素矩阵各变量表达式第44页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日455非定常不可压缩流动与自然对流换热问题求解不可压缩流动控制方程封闭空腔内自然对流换热基本控制方程解非定常不可压缩Navier-Stokes方程的时间分裂法时间分裂法求解自然对流问题时间分裂法求解Navier-Stokes方程的误差分析第45页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日465.1非定常不可压缩流动控制方程非定常不可压缩流动无因次Navier-Stokes方程为，以顶盖移动的方腔驱动流为例，定义无因次参量如下，第46页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日475.2封闭空腔内自然对流换热基本控制方程封闭空腔内自然对流换热基本控制方程第47页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日48无因次量定义对正方形封闭空腔，无因次参数的定义如下，对空气来说，第48页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日495.3解非定常不可压缩Navier-Stokes方程的时间分裂法Navier-Stokes可写成如下形式，对时间积分得，设为三步中间值，则分三步求解的过程为：非线性步、压力步、粘性步第49页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日50非线性步非线性步使用显式4阶Runge-Kutta方法，替代传统时间分裂法的Adams-Bashforth方法，使本时间步仅与上一时间步有关，此步中代入初始条件，不使用边界条件。第50页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日51压力步假设满足不可压缩条件，且在边界上法向上的投影为零，即，得,此步为解带有第二类边界条件的压力Poisson方程，解出压力Poisson方程后，再由上式求解。第51页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日52粘性步利用Crank-Nicolson格式离散粘性步，得Helmholtz方程，施以固体壁面边界条件即可求解第52页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日535.4时间分裂法求解自然对流问题求解自然对流问题的过程同Navier-Stokes方程，也是将一个时间步分为三步求解，所不同的是在一个时间步内除了解Navier-Stokes外，还需求解能量方程，能量方程中没有压力项，因此在第二步压力Poisson方程求解中与能量方程无关。第53页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日54非线性步非线性步仍使用显式4阶Runge-Kutta方法，第54页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日55压力步为解带有第二类边界条件的压力Poisson方程，解出压力Poisson方程后，再由(5.34)式求解。第55页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日56粘性步利用Crank-Nicolson格式离散粘性步，得Helmholtz方程，施以固体壁面边界条件即可求解。第56页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日575.5时间分裂法求解Navier-Stokes方程的误差分析利用时间分裂法求解Navier-Stokes方程基于以下两个假设，

满足不可压缩约束，即满足连续性方程，在边界的法线上的投影为零，高精度压力边界条件描述，压力Poisson方程变为，第57页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日586面向对象谱元方法程序设计面向对象程序设计思想面向对象谱元方法程序设计节点类标准单元类单元类整体区域类谱元方法程序设计中的几个问题插值函数积分的存储节点的编码及总体矩阵的形成Navier-Stokes方程和能量方程的求解过程第58页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日596.1面向对象程序设计思想按人们通常的思维方式建立模型，设计尽可能自然地表现求解方法的软件。为此，必须建立直接表现问题中事物与事物相互联系的概念，建立适合人们思维方式的描述方法；在面向对象的方法中主要有以下概念：对象、消息传递、类和继承、方法等。对象和消息传递是表示事物与事物间的关系；类和继承是适应人们思维方式对事物的描述方法；方法是作用在对象上的各种操作。第59页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日606.2面向对象谱元方法程序设计面向对象的程序设计包括面向对象分析、面向对象设计和面向对象实现对谱元方法来说，从谱元方法分析问题的思路出发，找出描述谱元方法的对象类（数据和方法的结合体），建立对象类的关系，并集成这些对象。设计节点类、标准单元、单元类、总体区域类。这些类中用属性来描述类的特性，用方法来实现某些操作。有些类中的属性可能又需要用对象类来描述，方法的作用对象也可能是某一类的实体。第60页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日61节点类主要的数据及方法包含节点的属性，如节点的编号、节点的位置坐标、节点的边界属性、微分方程右端项及边界条件值、所要求解的变量及其导数值。节点类定义，classCGridPoint{protected://节点的属性标准节点，编码，坐标，边界属性，方程右端项；压力，速度及其导数；public://节点类方法属性设置，边界设置，压力和速度设置；属性取值，边界属性取值；}第61页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日62标准单元类包括了在各求解单元中的常量。类中的数据及方法主要有：单元内网格的划分(含两方向的节点数及总节点数)、标准插值节点、插值函数的导数矩阵及有关标准插值函数的积分。标准单元类定义，classCStdElement{protected://标准单元类的属性网格属性，标准插值节点，插值函数的导数，插值函数积分；public://标准单元类方法标准单元类的构造，Chebyshev多项式及其导数；插值函数的积分，属性设置及取值；}第62页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日63单元类主要数据及方法包含元素序号、单元内节点划分(和标准单元内划分相对应)、标准单元类、节点类和总体坐标与局部坐标的变换关系。单元类定义，classCElement{protected://单元类的属性单元及网格属性，标准单元类，坐标变换，单元矩阵；public://单元类方法单元构造，标准单元类的构造，节点类构造；元素矩阵的形成，属性设置及取值；}第63页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日64总体区域类主要数据及方法包含区域内单元数、总节点数、节点类、标准单元类、总体矩阵等。总体区域类定义，classCDomain{protected://整体区域类的属性区域的属性，单元及网格属性，单元类，标准单元类；总体矩阵，对角元位置，边界节点序号；public://整体区域类方法单元构造，标准单元类的构造，节点类构造；元素矩阵的形成，属性设置及取值；}第64页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日656.2谱元方法程序设计中的几个问题插值函数积分的存储节点的编码及总体矩阵的形成对称带状矩阵的一维存储元素的编码影响元素和影响节点对角元的位置及一维数组的长度总刚度矩阵第一类边界条件的处理第65页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日66插值函数积分的存储在形成元素矩阵时，要计算如下积分，

可以看出，

在(i,j,k)变化时，不必占用(Nx+1)3个存储单元，这样既节约了存储空间又节省了计算时间。Bijk所占的空间可由下式计算，

Bijk的检索按如下方式进行，先将i,j,k按由大到小排列，设i>j>k，则Bijk在数组中的位置为，第66页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日676.3Navier-Stokes方程和能量方程的求解过程构造标准单元stdElement，求出对各单元不变的量。构造求解域pGeometry，对求解域进行网格划分，同时构造节点类pDot[NG]，构造各求解单元pElement[Ne]，对节点类中U、V、Φ赋初值，赋零值，Φ在热边界上赋值1，在冷边界上赋0，其余两边为第二类边界条件，在形成单元矩阵时自动满足。利用Runge-Kutta法求解、、，将求解结果存于、、解压力Poisson方程，第二类边界条件，解出压力后更新、、求解关于、、的Helmholtz方程，边界条件第一或第二类。解出Helmholtz方程后更新、、，即完成了一个时间步的求解。第67页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日687计算实例常微分方程求解Helmholtz和Poisson的求解简单边界Helmholtz方程求解复杂边界Poisson的求解移动顶盖方腔驱动非定常流动封闭方腔中的自然对流第68页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日697.1常微分方程求解【例1】用谱方法求解下面方程该问题的精确解的表达式为：

将求解域[0,-1]等分为10个求解单元,单元内使用5阶Chebyshev多项式进行谱逼近。数值计算结果与精确解的比较如下,比较表明谱元方法具有相当高的精度。第69页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日70例1数值计算结果与精确解的比较如下序号x数值解解析解绝对误差0123456789100.00.10.20.30.40.50.60.70.80.91-3.83295107661418e-0180.1559246312085910.3082929811551880.4535843017646150.5883485562994740.7092408909868030.8130550556913160.8967554397068940.957507402575330.99270559687679100.1559246312085960.3082929811551960.4535843017646250.5883485562994880.7092408909868240.8130550556913390.8967554397069170.9575074025753540.9927055968768061-3.83295107661418e-018-4.88498130835069e-015-8.04911692853239e-015-9.99200722162641e-015-1.45439216225896e-014-2.04281036531029e-014-2.33146835171283e-014-2.38697950294409e-014-2.34257058195908e-014-1.62092561595273e-014-2.22044604925031e-016第70页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日71【例2】【例2】用谱元方法求解下面方程该问题的精确解的表达式为：

将求解域[0,-1]等分为4个求解单元,单元内使用5阶Chebyshev多项式进行谱逼近。数值计算结果与精确解的比较如下,比较表明谱元方法具有相当高的精度第71页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日72例2数值计算结果与精确解的比较如下0123456789100.00.10.20.30.40.50.60.70.80.91-8.3294087758479e-0190.0147662965508840.0286795455946950.04087816195710.0504833997185180.0565905580149630.0582599254786710.0545073764403450.0442945201363430.0265183090547675.337544225295e-018

00.0147662965612980.0286795455705450.0408781618890680.050483399757920.0565905580149630.0582599255415590.0545073763178480.0442945199587630.0265183091540420

-8.3294087758479e-019-1.041408452529e-0112.4150008381962e-0116.8031476285757e-011-3.9402599238958e-011-3.816391647149e-016-6.2888715068876e-0111.2249663089436e-0101.7757908338245e-010-9.9275830611756e-0115.337544225295e-018

第72页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日737.2简单边界Helmholtz方程求解具体考察Helmholtz方程，精确解为，采用矩形单元采用曲边四边形单元第73页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日74两种单元求解结果比较(简单边界)表7-1(采用矩形单元的求解结果)表7-2(采用曲边四边形单元的求解结果)第74页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日75复杂边界Poisson的求解考察方程边界条件由精确解给出采用环形段单元采用曲边四边形单元第75页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日76两种单元求解结果比较(复杂边界)表7-3(采用环形段单元的求解结果)表7-4(采用曲边四边形单元的求解结果)第76页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日777.3移动顶盖方腔驱动非定常流动是计算流体力学典型数值实验模型，在很多情况下都是用它来检验数值方法的有效性；除移动顶盖外，初始速度场为零；移动顶盖上速度u=1,v=0，其余边界速度为零；取9个求解单元，使用8阶Chebyshev多项式插值；计算了Re=100的情况，取两种时间步长T1=0.004，T2=0.005，计算2000和1600步，各需12和9.6小时第77页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日78移动顶盖方腔驱动流物理模型第78页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日79竖直中心线上的速度U,水平中心线上的速度V图7-3竖直中心线上的速度U图7-4水平中心线上的速度V第79页，共94页，2022年，5月20日，8点50分，星期三2022年9月13日80方腔内中心点的速度