谱元方法课程ppt课件.ppt

《谱元方法课程ppt课件.ppt》由会员分享，可在线阅读，更多相关《谱元方法课程ppt课件.ppt（94页珍藏版）》请在三一办公上搜索。

1、谱元方法简介,主讲：秦国良电话：82663537邮箱：,2022年11月16日,Tel：82663537 e-mail：,2,1、正交函数系与谱近似,1.1 正交函数系与正交多项式1.2 函数的Fourier展开1.3 Chebyshev谱逼近(离散函数的Fourier展开)1.4 插值函数的导数1.5 二维函数的Chebyshev谱近似,2022年11月16日,Tel：82663537 e-mail：,3,1.1 正交函数系与正交多项式,设函数系j(x)，如满足条件， 则称函数系j(x)，在区间a,b上关于权函数w(x) 正交，当j(x)是代数多项式时，称为正交多项式。如Chebyshev多

2、项式，,2022年11月16日,Tel：82663537 e-mail：,4,Chebyshev多项式的正交性,2022年11月16日,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展开定义为，,2022

3、年11月16日,Tel：82663537 e-mail：,6,1.3 Chebyshev谱逼近(离散函数的Fourier展开),取节点系xk, k = 0,1,N为N阶Chebyshev多项式的极值点(称Gauss-Lobatto点)，即，权系数则上面的展开式变为，此式亦可看作f(x)在配置点上的插值, g(x)为插值函数。,2022年11月16日,Tel：82663537 e-mail：,7,1.4 插值函数的导数,对 微分得：插值函数在配置点的导数可表示为，,2022年11月16日,Tel：82663537 e-mail：,8,1.5 二维函数的Chebyshev谱近似,对二维函数u(,)

4、，在标准正方形单元内，定义,方向节点系j，j=0,1,Nx和k ，k=0,1,Ny，,2022年11月16日,Tel：82663537 e-mail：,9,2 谱方法求解微分方程,2.1 解微分方程的加权残量法MWR （Method of Weight Residuals)2.2 谱方法求解微分方程 2.3 小结,2022年11月16日,Tel：82663537 e-mail：,10,2.1 解微分方程的加权残量法MWR（Method of Weight Residuals),基本思想MWR的步骤MWR的基本方法,2022年11月16日,Tel：82663537 e-mail：,11,加权残量

5、法基本思想(1),设微分方程边值问题为，加权余量法，是求微分方程形式如下的近似解。a=a0,a1,aNT为待定向量，0, 1, N为为基函数，且是一组线性无关的函数,2022年11月16日,Tel：82663537 e-mail：,12,加权残量法基本思想(2),“残量”定义如下，显然只有当试函数为边值问题的精确解时，余量(2.31)才为零。加权残量法适当选取待定向量a=a0,a1,aNT，使得“残量”极小。通常某一权函数系wi(x),i=0,1,N，使得权函数与“残量”正交，来确定待定向量。,2022年11月16日,Tel：82663537 e-mail：,13,MWR的步骤,选取一组满足要

6、求的基函数构造试函数 并满足所有边界条件；选取一组权函数运用MWR准则， 得到 的代数方程组； 解上述代数方程组，确定待定参数。,2022年11月16日,Tel：82663537 e-mail：,14,MWR的基本方法,按权函数的不同有以下几种基本方法Galerkin法(Galerkin Method)配置法(Collocation Method)Tan方法(Tan Method) 类似于Galerkin法，Tan方法中，试函数不需要满足边界条件，增加边界约束系数来满足边界条件,2022年11月16日,Tel：82663537 e-mail：,15,2.2 谱方法求解微分方程,以非线性Burg

7、er方程为例来说明谱方法的离散过程,并区别Galerkin法,Tau方法和Collocation方法。Burger方程如下:,2022年11月16日,Tel：82663537 e-mail：,16,Chebyshev Galerkin方法 (Galerkin Method,方程的近似解表达为， 由MWR准则得,2022年11月16日,Tel：82663537 e-mail：,17,Chebyshev Tau方法(Tan Method,方程的近似解表达为， 由MWR准则得, 边界约束,2022年11月16日,Tel：82663537 e-mail：,18,Chebyshev Collocatio

8、n方法(Collocation Method),选取节点系为N阶Chebyshev多项式的极值点, 方程的近似解用下式表示，由Collocation方程得,2022年11月16日,Tel：82663537 e-mail：,19,2.3 小结,一种求解偏微分方程的高阶精度数值方法。 属于求解偏微分方程加权残量法的特例。谱方法使用高阶正交多项式作为展开函数 。谱方法最受人青睐的优越性在于它具有“无穷阶”的收敛速度，其确切含义为，若原微分方程的解无限可微，则由适当的谱方法所得到的近似解对原问题的收敛速度比1/N的任何幂都更快，这里N是所取基函数的个数。Kreiss和Oliger, OrszagGot

9、tlieb和Orszag,2022年11月16日,Tel：82663537 e-mail：,20,3 谱元方法求解微分方程,回顾谱方法 使用加权残量法谱元方法基本思想常微分方程边值问题的Galerkin变分原理偏微分方程的Galerkin变分原理Galerkin逼近解常微分方程元素矩阵的形成偏微分方程元素矩阵的形成极坐标系下环形区域的谱元方法总刚度矩阵的形成,2022年11月16日,Tel：82663537 e-mail：,21,3.1 回顾谱方法 使用加权残量法,选取一组满足要求的基函数,构造试函数，满足所有的边界条件,选取一组权函数,利用MWR准则，得到ai 的代数方程,求解上述代数方程组

10、，确定待定参数,2022年11月16日,Tel：82663537 e-mail：,22,3.2 谱元方法基本思想,使用有限元方法的思想 ，将求解域分成若干子域 采用谱方法中谱近似技术代替有限元中的插值函数 采用有限元中等参单元,亦可逼近复杂求解边界,2022年11月16日,Tel：82663537 e-mail：,23,3.3 常微分方程边值问题的Galerkin变分原理,考察二点边值问题上述问题的Galerkin变分问题是：求 ，使得，其中,2022年11月16日,Tel：82663537 e-mail：,24,3.4 偏微分方程的Galerkin变分原理,考察二维椭圆边值问题 上述问题的G

11、alerkin变分问题是：求 ，使得，,2022年11月16日,Tel：82663537 e-mail：,25,3.5 Galerkin逼近解,设Vh是V的有限维子空间，当h-0时，Vh的维数无限增加，直到充满V为止，那么Galerkin逼近解 ，使得，设 是Vh的基函数系，那么由于 ，故有， 其中得由于 的任意性，得，即得线性代数方程组，可以证明此代数方程组的解唯一存在，即可解出Galerkin逼近解。,2022年11月16日,Tel：82663537 e-mail：,26,3.6 常微分方程元素矩阵的形成,考察以下微分方程用谱方法求解此微分方程,先将求解区域分成若干单元,设单元总数为N,考

12、察第i个单元ei,将ei转化为标准求解单元,设x为总体坐标,为标准求解单元内的局部坐标,并设,2022年11月16日,Tel：82663537 e-mail：,27,插值点和插值函数,在第i个元素内选取Nx+1个点作为u(x)和f(x)的插值点,本文选取阶Chebyshev多项式的极值点作为插值点,即 ,u(x)和f(x)可以表达为，,2022年11月16日,Tel：82663537 e-mail：,28,元素矩阵形成,对一维边值问题,在ei单元内线性方程式变为 其中，上面线性方程组变为，,2022年11月16日,Tel：82663537 e-mail：,29,的推导过程,2022年11月16

13、日,Tel：82663537 e-mail：,30,3.7 偏微分方程元素矩阵的形成,具体考察Helmholtz方程，其中，作以下坐标变换在标准正方形单元内，使用伪谱Chebyshev谱逼近u(,)，则元素内线性方程组式变为，,2022年11月16日,Tel：82663537 e-mail：,31,总体坐标系与局部坐标系,2022年11月16日,Tel：82663537 e-mail：,32,元素矩阵,2022年11月16日,Tel：82663537 e-mail：,33,3.8 极坐标系下环形区域的谱元方法,在极坐标系下，Helmholtz方程表达为，将极坐标系下环形段单元，变换为直角坐标系

14、下的正方形单元，坐标变换关系如下，元素矩阵的形式为，,2022年11月16日,Tel：82663537 e-mail：,34,极坐标坐标系与局部坐标系,2022年11月16日,Tel：82663537 e-mail：,35,元素矩阵,式中各变量的表达式如下，,2022年11月16日,Tel：82663537 e-mail：,36,3.9 总刚度矩阵的形成,节点影响元素集，节点所在的元素的并集影响节点，对某点，其影响元素内的所有节点相互影响元素集，两个节点和影响元素之交集对节点和，总体矩阵元素为，由于 在一个元素e内计算，故称它为单元矩阵，上式表明矩阵有贡献的只是()影响元素集所有单元矩阵相应的

15、元素。,2022年11月16日,Tel：82663537 e-mail：,37,4 复杂区域的谱元方法,采用谱逼近的等参数单元总体坐标系下导数的计算坐标变换矩阵及变换行列式局部坐标系下偏导数的计算总体坐标系下偏导数的计算复杂求解区域的谱元方法,2022年11月16日,Tel：82663537 e-mail：,38,总体坐标与局部坐标系,2022年11月16日,Tel：82663537 e-mail：,39,4.1采用谱逼近的等参数单元,总体坐标与局部坐标(如图所示)等参单元，插值函数表达式和坐标变换表达式一致函数的插值表达为，坐标变换的表达式为，导数表达式,2022年11月16日,Tel：82

16、663537 e-mail：,40,4.2 总体坐标系下导数的计算,坐标变换矩阵(Jacobi矩阵)及变换行列式(Jacobi行列式),微元面积、微元弧长 (为常数)时， (为常数)时,2022年11月16日,Tel：82663537 e-mail：,41,局部坐标系下偏导数的计算,由插值式对求偏导得，上式在(m,n)处的值为，比照函数得，则局部坐标系下的偏导数表达为，,2022年11月16日,Tel：82663537 e-mail：,42,总体坐标系下偏导数的计算,由插值函数得，由Jacobi矩阵得总体坐标系下，偏导数计算式为，,2022年11月16日,Tel：82663537 e-mail

17、：,43,4.3 复杂求解区域的谱元方法,考察Helmholtz方程则元素矩阵为，,2022年11月16日,Tel：82663537 e-mail：,44,元素矩阵各变量表达式,2022年11月16日,Tel：82663537 e-mail：,45,5 非定常不可压缩流动与自然对流换热问题求解,不可压缩流动控制方程封闭空腔内自然对流换热基本控制方程解非定常不可压缩Navier-Stokes方程的时间分裂法时间分裂法求解自然对流问题时间分裂法求解Navier-Stokes方程的误差分析,2022年11月16日,Tel：82663537 e-mail：,46,5.1 非定常不可压缩流动控制方程,非

18、定常不可压缩流动无因次Navier-Stokes方程为，以顶盖移动的方腔驱动流为例，定义无因次参量如下，,2022年11月16日,Tel：82663537 e-mail：,47,5.2 封闭空腔内自然对流换热基本控制方程,封闭空腔内自然对流换热基本控制方程,2022年11月16日,Tel：82663537 e-mail：,48,无因次量定义,对正方形封闭空腔，无因次参数的定义如下，对空气来说，,2022年11月16日,Tel：82663537 e-mail：,49,5.3 解非定常不可压缩Navier-Stokes方程的时间分裂法,Navier-Stokes可写成如下形式，对时间积分得，设 为

19、三步中间值，则分三步求解的过程为：非线性步、压力步、粘性步,2022年11月16日,Tel：82663537 e-mail：,50,非线性步,非线性步使用显式4阶Runge-Kutta方法，替代传统时间分裂法的Adams-Bashforth方法，使本时间步仅与上一时间步有关，此步中代入初始条件，不使用边界条件。,2022年11月16日,Tel：82663537 e-mail：,51,压力步,假设 满足不可压缩条件，且在边界上法向上的投影为零，即，得,此步为解带有第二类边界条件的压力Poisson方程，解出压力Poisson方程后，再由上式求解 。,2022年11月16日,Tel：8266353

20、7 e-mail：,52,粘性步,利用Crank-Nicolson格式离散粘性步，得Helmholtz方程，施以固体壁面边界条件即可求解,2022年11月16日,Tel：82663537 e-mail：,53,5.4 时间分裂法求解自然对流问题,求解自然对流问题的过程同Navier-Stokes方程，也是将一个时间步分为三步求解，所不同的是在一个时间步内除了解Navier-Stokes外，还需求解能量方程，能量方程中没有压力项，因此在第二步压力Poisson方程求解中与能量方程无关。,2022年11月16日,Tel：82663537 e-mail：,54,非线性步,非线性步仍使用显式4阶Run

21、ge-Kutta方法，,2022年11月16日,Tel：82663537 e-mail：,55,压力步,为解带有第二类边界条件的压力Poisson方程，解出压力Poisson方程后，再由(5.34)式求解 。,2022年11月16日,Tel：82663537 e-mail：,56,粘性步,利用Crank-Nicolson格式离散粘性步，得Helmholtz方程，施以固体壁面边界条件即可求解。,2022年11月16日,Tel：82663537 e-mail：,57,5.5 时间分裂法求解Navier-Stokes方程的误差分析,利用时间分裂法求解Navier-Stokes方程基于以下两个假设，

22、满足不可压缩约束，即满足连续性方程， 在边界的法线上的投影为零，高精度压力边界条件描述，压力Poisson方程变为，,2022年11月16日,Tel：82663537 e-mail：,58,6 面向对象谱元方法程序设计,面向对象程序设计思想面向对象谱元方法程序设计节点类标准单元类单元类整体区域类谱元方法程序设计中的几个问题插值函数积分的存储节点的编码及总体矩阵的形成Navier-Stokes方程和能量方程的求解过程,2022年11月16日,Tel：82663537 e-mail：,59,6.1 面向对象程序设计思想,按人们通常的思维方式建立模型，设计尽可能自然地表现求解方法的软件。为此，必须建

23、立直接表现问题中事物与事物相互联系的概念，建立适合人们思维方式的描述方法；在面向对象的方法中主要有以下概念：对象、消息传递、类和继承、方法等。对象和消息传递是表示事物与事物间的关系；类和继承是适应人们思维方式对事物的描述方法；方法是作用在对象上的各种操作。,2022年11月16日,Tel：82663537 e-mail：,60,6.2 面向对象谱元方法程序设计,面向对象的程序设计包括面向对象分析、面向对象设计和面向对象实现对谱元方法来说，从谱元方法分析问题的思路出发，找出描述谱元方法的对象类（数据和方法的结合体），建立对象类的关系，并集成这些对象。设计节点类、标准单元、单元类、总体区域类。这些

24、类中用属性来描述类的特性，用方法来实现某些操作。有些类中的属性可能又需要用对象类来描述，方法的作用对象也可能是某一类的实体。,2022年11月16日,Tel：82663537 e-mail：,61,节点类,主要的数据及方法包含节点的属性，如节点的编号、节点的位置坐标、节点的边界属性、微分方程右端项及边界条件值、所要求解的变量及其导数值。节点类定义，class CGridPoint protected:/节点的属性 标准节点，编码，坐标，边界属性，方程右端项； 压力，速度及其导数；public:/节点类方法 属性设置，边界设置，压力和速度设置； 属性取值，边界属性取值；,2022年11月16日,

25、Tel：82663537 e-mail：,62,标准单元类,包括了在各求解单元中的常量。类中的数据及方法主要有：单元内网格的划分(含两方向的节点数及总节点数)、标准插值节点、插值函数的导数矩阵及有关标准插值函数的积分。标准单元类定义，class CStdElement protected:/标准单元类的属性 网格属性，标准插值节点，插值函数的导数，插值函数积分；public:/标准单元类方法 标准单元类的构造，Chebyshev多项式及其导数； 插值函数的积分，属性设置及取值；,2022年11月16日,Tel：82663537 e-mail：,63,单元类,主要数据及方法包含元素序号、单元内节

26、点划分(和标准单元内划分相对应)、标准单元类、节点类和总体坐标与局部坐标的变换关系。单元类定义，class CElement protected:/单元类的属性 单元及网格属性，标准单元类，坐标变换，单元矩阵；public:/单元类方法 单元构造，标准单元类的构造，节点类构造； 元素矩阵的形成，属性设置及取值；,2022年11月16日,Tel：82663537 e-mail：,64,总体区域类,主要数据及方法包含区域内单元数、总节点数、节点类、标准单元类、总体矩阵等。总体区域类定义，class CDomain protected:/整体区域类的属性 区域的属性，单元及网格属性，单元类，标准单元

27、类； 总体矩阵，对角元位置，边界节点序号；public:/整体区域类方法 单元构造，标准单元类的构造，节点类构造； 元素矩阵的形成，属性设置及取值；,2022年11月16日,Tel：82663537 e-mail：,65,6.2 谱元方法程序设计中的几个问题,插值函数积分的存储节点的编码及总体矩阵的形成对称带状矩阵的一维存储元素的编码影响元素和影响节点对角元的位置及一维数组的长度总刚度矩阵第一类边界条件的处理,2022年11月16日,Tel：82663537 e-mail：,66,插值函数积分的存储,在形成元素矩阵时，要计算如下积分， 可以看出， 在(i,j,k)变化时，不必占用(Nx+1)3

28、个存储单元，这样既节约了存储空间又节省了计算时间。Bijk所占的空间可由下式计算， Bijk的检索按如下方式进行，先将i,j,k按由大到小排列，设ijk，则Bijk在数组中的位置为，,2022年11月16日,Tel：82663537 e-mail：,67,6.3 Navier-Stokes方程和能量方程的求解过程,构造标准单元stdElement，求出对各单元不变的量。构造求解域pGeometry，对求解域进行网格划分，同时构造节点类pDotNG，构造各求解单元pElementNe，对节点类中U、V、赋初值，赋零值，在热边界上赋值1，在冷边界上赋0，其余两边为第二类边界条件，在形成单元矩阵时自

29、动满足。利用Runge-Kutta法求解 、 、 ，将求解结果存于 、 、 解压力Poisson方程，第二类边界条件，解出压力后更新 、 、求解关于 、 、 的Helmholtz方程，边界条件第一或第二类。解出Helmholtz方程后更新 、 、 ，即完成了一个时间步的求解。,2022年11月16日,Tel：82663537 e-mail：,68,7 计算实例,常微分方程求解Helmholtz和Poisson的求解简单边界Helmholtz方程求解复杂边界Poisson的求解移动顶盖方腔驱动非定常流动封闭方腔中的自然对流,2022年11月16日,Tel：82663537 e-mail：,69,

30、7.1 常微分方程求解,【例1】用谱方法求解下面方程该问题的精确解的表达式为： 将求解域0,-1等分为10个求解单元,单元内使用5阶Chebyshev多项式进行谱逼近。数值计算结果与精确解的比较如下,比较表明谱元方法具有相当高的精度。,2022年11月16日,Tel：82663537 e-mail：,70,例1 数值计算结果与精确解的比较如下,2022年11月16日,Tel：82663537 e-mail：,71,【例2】,【例2】用谱元方法求解下面方程该问题的精确解的表达式为： 将求解域0,-1等分为4个求解单元,单元内使用5阶Chebyshev多项式进行谱逼近。数值计算结果与精确解的比较如

31、下,比较表明谱元方法具有相当高的精度,2022年11月16日,Tel：82663537 e-mail：,72,例2 数值计算结果与精确解的比较如下,2022年11月16日,Tel：82663537 e-mail：,73,7.2 简单边界Helmholtz方程求解,具体考察Helmholtz方程，精确解为，采用矩形单元采用曲边四边形单元,2022年11月16日,Tel：82663537 e-mail：,74,两种单元求解结果比较(简单边界),表7-1(采用矩形单元的求解结果)表7-2(采用曲边四边形单元的求解结果),2022年11月16日,Tel：82663537 e-mail：,75,复杂边界

32、Poisson的求解,考察方程边界条件由精确解给出采用环形段单元采用曲边四边形单元,2022年11月16日,Tel：82663537 e-mail：,76,两种单元求解结果比较(复杂边界),表7-3(采用环形段单元的求解结果)表7-4(采用曲边四边形单元的求解结果),2022年11月16日,Tel：82663537 e-mail：,77,7.3 移动顶盖方腔驱动非定常流动,是计算流体力学典型数值实验模型，在很多情况下都是用它来检验数值方法的有效性；除移动顶盖外，初始速度场为零；移动顶盖上速度u=1,v=0，其余边界速度为零；取9个求解单元，使用8阶Chebyshev多项式插值；计算了Re=10

33、0的情况，取两种时间步长T1=0.004，T2=0.005，计算2000和1600步，各需12和9.6小时,2022年11月16日,Tel：82663537 e-mail：,78,移动顶盖方腔驱动流物理模型,2022年11月16日,Tel：82663537 e-mail：,79,竖直中心线上的速度U, 水平中心线上的速度V,图7-3 竖直中心线上的速度U,图7-4 水平中心线上的速度V,2022年11月16日,Tel：82663537 e-mail：,80,方腔内中心点的速度U 和 速度V随时间的变化,图7-5 方腔内中心点的速度U,图7-6 方腔内中心点的速度V,2022年11月16日,Te

34、l：82663537 e-mail：,81,封闭方腔中的自然对流,计算中网格单元数为3x3，网格内部的插值函数采用8阶Chebyshev多项式；计算的空气在正方形封闭空腔内的自然对流，空气的Pr数取为0.71Ra=103, 时间步长取T=3x10-4，计算1500时间步；耗时12小时Ra=104 , 时间步长取T=3x10-4，计算1500时间步；耗时12小时Ra=105 , 时间步长取T=4x10-5，计算12500时间步；耗时90小时,2022年11月16日,Tel：82663537 e-mail：,82,封闭方腔中的自然对流物理模型,2022年11月16日,Tel：82663537 e-

35、mail：,83,数值模拟结果与基准解的比较,数值模拟结果与基准解的比较,2022年11月16日,Tel：82663537 e-mail：,84,竖直中心线u速度分布、水平中心线上v速度分布,图7-13 竖直中心线u速度分布,图7-14水平中心线上v速度分布,2022年11月16日,Tel：82663537 e-mail：,85,水平中心线上的温度分布,图7-15水平中心线上的温度分布,2022年11月16日,Tel：82663537 e-mail：,86,Ra=103速度矢量、温度场,2022年11月16日,Tel：82663537 e-mail：,87,Ra=104的速度矢量和温度场,20

36、22年11月16日,Tel：82663537 e-mail：,88,Ra=105的速度矢量和温度场,总结,函数逼近谱元方法和有限元方法时间分裂法解N-S方程发展谱元方法的几个问题,2022年11月16日,Tel：82663537 e-mail：,89,2022年11月16日,Tel：82663537 e-mail：,90,了解谱元方法,谱元,谱方法,有限元方法,Spectral,Element,2022年11月16日,Tel：82663537 e-mail：,91,有限元类方法求解过程,方程变分,区域划分,定基函数,单元分析,合成总刚,边界条件,问题的解,原始方程,2022年11月16日,Tel：82663537 e-mail：,92,谱元方法的基本原理,2022年11月16日,Tel：82663537 e-mail：,93,时间分裂法的应用,2022年11月16日,Tel：82663537 e-mail：,94,谱元方法的若干问题,时间方向的精度问题（高Re数问题）,复杂区域的进一步适应性,计算速度问题,圆形区域问题,一般区域问题,极坐标或者柱坐标系,等参谱元方法,提高时间分裂法的阶数,提高压力边界条件阶数,