计算流体力学(中科院力学所)-第4讲-差分方法2

计算流体力学讲义第四讲有限差分法（2）李新亮[e-mail]；力学所主楼219；82543801知识点：离散误差的Fourier分析;间断周围数值振荡的原因；GVC格式；模型方程向N-S方程的推广；1讲义、课件上传至(流体中文网）-“流体论坛”-“CFD基础理论”CopyrightbyLiXinliang§3.3差分格式的进一步分析1.耗散与色散误差2CopyrightbyLiXinliang)sin()0,(0xxuxutu精确解1阶迎风xu0123456-1-0.500.51Exact1stupwind2ndupwindCopyrightbyLiXinliangt=solutionsof1DconvectionEq.2阶迎风)2/()43(21xuuuujjjjxuuujjj/)(1数值实验boundaryperiodicxxxuxutu],2,0[)sin()0,(0时间推进：3步TVD型Runge-Kutta,且时间步长足够小（误差忽略）空间离散：1阶及2阶迎风格式（20个网格点）实验观察到的现象——两类误差：振幅误差相位误差（波速误差）CopyrightbyLiXinliang3对以上“实验现象”进行理论分析半离散分析：假设时间推进是精确的，仅分析空间离散带来的误差（难度小、常用）全离散分析：同时分析时、空离散的误差（难度大）boundaryperiodicxexfcxfctfikx],2,0[)0,()0(0考查问题：实际上就是普通三角函数，采用复数形式仅仅是为了理论推导方便。用实数形式sin(kx),cos(kx)推导形式上略显繁琐。精确解：ikxikctctxikeeetxf)(),(差分格式：0jxjfctf（1）xfffjjjx/)(1)2/()43(21xuuufjjjjx其他格式……jikxjexu)(假设对于：有jikxjxexku~隐含假设:线性差分格式非线性系统作用于单波，会产生多个谐波（2）差分没有误差xikk~CopyrightbyLiXinliang4令：jikxjetutxf)(ˆ),(jikxjxxektutxf/~)(ˆ),(（1）式化为：0/~)(ˆ)(ˆjjikxikxxektucedttud0/~)(ˆ)(ˆxktucdttud“半离散化”：空间导数差分计算，时间方程（常微）精确计算xctkeutu/~)0(ˆ)(ˆxctkikxikxjjjeuetutxf/~)0(ˆ)(ˆ),(ixikk~如果，无误差分析（修正波数）与误差的关系k~irikkk~))/((/)0(ˆ),(xkctkxikxctkjijreeutxf理想情况：的误差导致解的幅值误差——耗散误差的误差导致解传播速度的误差——色散误差xkkkir,0rkikjikxjexu)(假设对于：有jikxjxexku~反映了一个波内的点数。PPW（波内的点数）=CopyrightbyLiXinliang5jikxjexu)(jikxjxexku~耗散、色散误差分别由修正波数的实部和虚部决定。k~关键参数：修正波数含义：反应波数（谱）空间内差分的误差jikxkkjefxfˆ)(任意函数：jikxkkjefikxfˆ)(定义：kkfikfˆˆ求导数，精确解差分解jikxkkjxjefxkifFˆ/~kkfxkFˆ~ˆFourier分析的任务计算出，并考差其与的逼近程度。k~ixik考察格式分辨率（resolution）的重要指标精度：反映时的情况分辨率：网格点数很少（例如波里面只有6个点）时的性能对于多尺度问题，分辨率更重要。牺牲精度，提高分辨率0xxk优秀的差分格式，1个波长里面6个点即可/2精度分辨率CopyrightbyLiXinliang6如何计算修正波数？jikxjexu)(jikxjxexku~定义：方法1.理论计算根据差分具体表达式及定义计算例1：xuuujjjx1令jikxjeu则：jjjjikxiikxxxikikxjjjxexkexeeexxuuu~)1()(1)(1xk于是：sin)cos1(1~ieki1阶迎风sin,cos1irkk例2：)2/()43(21xuuuujjjjx2阶迎风)43(2243221iiikxjjjjxeexexuuuuj2/)43(~2iieek2/)2sinsin4(,2/)2coscos43(irkkCopyrightbyLiXinliang7方法2：数值计算jikxjexu)(jikxjxexku~定义：Step1）选取计算域[0,2],计算网格（例如64，128）Step2）给定波数k,生成函数值Step3)调用差分子程序，得到导数值Step4)通过Fourier反变换，得到谱：假设已有求差分的子程序（黑箱，已知是线性的）x线性黑箱ju强调：研究CFD本身，不能只使用理论手段，还要用数值手段jxjuv)...2,1(Njeujikxj)...2,1(NjvjNjikxjkjevNv11ˆ根据修正波数的定义，有xvkkˆ~Step5)改变k的值，重复2-5，得到对于的依赖关系。画图k~xk非线性情况会产生高次谐波，造成step4中隐含的假设无法成立将Fourier分析手段拓展到非线性系统需要研究的课题隐含条件：只有波数为k的那个谱不为0（线性系统）CopyrightbyLiXinliang8中心差分格式的色散特性0：精确解；1：4阶普通2：6阶普通；3：4阶紧致4：6阶紧致；5：6阶超紧致迎风差分格式的色散特性0：精确解，1：2阶迎风2：5阶迎风偏心3：3阶迎风紧致4：5阶迎风紧致每个波长里面2个网格点，谱方法的分辨率，差分法分辨率的极限（只有无穷阶精度才能达到）2)int(WavePersPoPPW20阶超紧致格式——接近谱方法CopyrightbyLiXinliang9不同差分格式的色散误差曲线结论：要求分辨率相同的情况下，采用高阶格式可放宽空间网格步长，从而减少计算量重要方向：高分辨率差分格式0：精确解1：2阶迎风2:3阶迎风3：3阶迎风紧致4：5阶迎风紧致指定误差要求的情况下，不同差分格式能模拟的最大(越大，所需网格越少）/2PPW作业题1：构造高分辨率差分格式，并进行理论分析及数值实验针对单波方程：0xutu对于空间导数，构造出一种不超过6点格式；并进行Fourier误差分析，画出kr,ki的曲线。要求：精度不限；网格基架点数不超过6个；能够分辨的波数范围尽量宽；（即kr,ki曲线近可能接近准确解）给出差分的具体表达式，画出kr,ki的曲线；说明构造格式的阶数，并采用本PPT第5页的方法给出的精度验证；26154131231jjjjjjjjuauauauauauaxuu16514233241jjjjjjjjuauauauauauaxuu形如：……另外，进行如下数值验证：)sin()0,(]2,0[,0xxuxxutu空间采用20个网格点，采用新构造的差分格式离散；时间推进采用3步Runge-Kutta方法，时间步长可足够小（例如0.01）。给出t=20,50两个时刻的数值解，与精确解比较（画图），并给出数值解的L2模误差。10CopyrightbyLiXinliang提示：1.如不使用优化技术，则格式构造方法简单，Taylor展开后解代数方程组即可。2.建议尝试使用优化技术26154131231jjjjjjjjuauauauauauaxuu例：假设格式形式如下如果要求其有5阶精度，则通过Taylor展开可得到6个方程，6个系数可直接解出。我们要求其有4阶精度（当然3阶，2阶也可），于是Taylor展开只能提供5个方程。6个未知数（a1-a6），5个方程；有1个自由参数。调整这个自由参数，使得kr,ki曲线最为理想。如何调整？1)可以人工调整，观察kr,ki曲线，选取满意的。2）可自动调整，设立一个优化目标函数。例如调整自由参数，使得该目标函数取最大值。思路：牺牲精度，提高分辨率05.0)(:***ik11CopyrightbyLiXinliang附录：部分差分格式…j-2j-1jj+1…0xuatu表中的迎风差分格式均针对a0当a0时，需把下标的“j+k”换成“j-k”（例如把j+2换成j-2,把j-1换成j+1);并在表达式前加上“-”号。例：迎风偏斜格式：上游的基架点更多些（或上游权重更大）xfffffjjjjj6/)236(112xfffffjjjjj6/)236(11212CopyrightbyLiXinliang§3.4数值解的群速度及间断处数值振荡来源jikxjexu)(对于：有jikxjxexku~修正波数)/(/)0(ˆ),(ijrctkxikxctkjeeutxfboundaryperiodicxexfcxfctfikx],2,0[)0,()0(0数值解色散误差：数值解传播的速度与精确解不一致1/ddki1/ddki数值解传播偏快数值解传播偏慢0：精确解；1：2阶迎风；2：5阶迎风偏心3：3阶迎风紧致；4：5阶迎风紧致快格式(FST)：慢格式(SLW):混合格式(MXD):1/ddki1/ddki特点：波数越高，误差越严重1.色散误差与群速度13CopyrightbyLiXinliang)sin(][)sin()0,(1)5.1(16)5.1(16016222xaeexaexuxxx0xutut=0.5时刻的精确解及数值解25,19410aa200/1xxuuffjjjj/)(2/)(11空间离散：五种不同格式；时间推进：3阶Runge-Kutta【数值实验】波的传播问题)2/()(11xuuujjj)2/()43('11xuuuujjjj)4/()('22xuuujjj观察现象：1）高波数成分误差严重，低波数成分误差不明显；2）二阶Pade格式的解传播速度快于精确解，其余格式偏慢；3）迎风型格式有耗散，尤其是二阶迎风格式；概念：群速度——波包传播的速度14CopyrightbyLiXinliang2.间断附近数值振荡的来源【数值实验】间断的传播1,0axuatu5.015.00)0,(xxxu计算域[0,1];计算网格点100时间推进：3阶Runge-Kutta空间离散：1）二阶中心差分)2/()43(21xuuuujjjjx)2/()(11xuuujjjx2/)2sinsin4(iksinikki00.511.522.5300.511.522ndcentrialscheme2ndupwindschemeki=2阶迎风及2阶中心格式的色散特性2）二阶迎风差分15CopyrightbyLiXinliang过激波数值振荡的根源——色散误差导致群速度不一致快格式慢格式波前振荡波后振荡=+++…群速度控制的基本思路（群速度控制GVC:Fu&Ma）：间断前、后分别采用快格式和慢格式，可有效抑制振荡Zhuang&Zhang:抑制波动原则示意图：间断的Fourier分解