关于期末考试

1、关于期末考试谢谢各位同学来听这门课！正如我所说，我们的期末考试是以大家做报告的形式进行。具体安排如下：计算天文期末汇报暨天文系第四届数值模拟研讨会时间：2009年1月9日下午2：006：00地点：天文系202会议室报告形式：口头（均利用ppt文件），每人（组）1015分钟联系人：P. F. Chen联系电话：83594651电子邮件：注意事项：注意事项：1.请大家好好准备，我会尽量邀请一些导师参加；2.可以单独做个题目，也可2人为一组合作做题目；3.题目可以是我提供的，也可结合你们各自的研究自行选择；4.除作报告外，每人（组）请将模拟工作整理出13页纸（形式如下），以便于我编辑成文集后发送给与

2、会人士：太阳耀斑环动力学数值模拟张三（南京大学天文系，南京 210093）1 引言他人工作等.2 数值方法.3 结果含图4 讨论/结论 参考文献注意事项：注意事项：5.ppt报告中请多用些图甚至动画来显示你的结果；6.期间有任何问题请及时与我(科技馆2415室)联系，讨论；7.请记住几个时间表： 1月1日前：请发信件将题目告诉我：我的电子邮件地址是 ; BBS帐号是MHD；两者任一。 1月8日前：所有同学将论文(pdf)及报告(ppt文件)拷贝到我办公室 （或通过邮件附件发给我）。 1月9日下午：天文系202室研讨会19. 高级算法之CIP法Yabe (矢部) et al. 1991, Com

3、put. Phys. Commun., 66, 219 2001, JCP, 169, 556发展了CIP (Cubic Interpolated Profile) 高精度差分法，举例如下：若 u为常数，则其解为 ，其差分形式为0 xfutf)0 ,(),(utxftxf),(),(ttuxfttxfii则新时刻的分布化为老时刻某空间点的分布ttt tu1ii1i其求解之最简方法为一次迎风格式（以u 0为例）)(11ninininiffxtuff但其缺点是强数值耗散。CIP方法在两格点间 采用3次多项式拟合),(1iixxniniiiiiiffxxxxbxxaxF)()()()(23拟合系数满

4、足此函数及其导数值在 i-1, i 格点处连续的条件:niifxF1)(1nixifxF1)(1的导数点处格表示为 fifni31212xffxffaninininiixffxffbninininii12123若u0 则 i-1 i+1， x - x 因此，CIP法要求同时求解为此，对原方程求导同样， 也满足综上所述，nif0 xfutfnif),(),(ttuxfttxfiininiiiiniffbatuxFf)(231niiiinifbatuxFxf23)(dd21tu初条由 和 同时给出 。nifnif计算实例：单波问题非线性方程1. 非对流项即上两式求解分两步实现Gxufgxfutfg

5、xuftfxufGxfutfGtfxufGtf2. 对流项0 xfutf0 xfutf1. 非对流项GtfxufGtf采用时间前差，空间中心差分tGffinii*txuufxffffftxuuftxGGffiinniiniiniiniiniiiii22)()(22111*11*11111*常可用交错网格1ixix1ix21ix21ixepepepuu2. 对流项采用CIP方法0 xfutf0 xfutf*231iiiiniffbaf*2123)(ddiiiinifbatuxFxf3*1*2*1*2xffxffaiiiiixffxffbiiiii*1*2*123tu若u0 则 i-1 i+1，

6、x - x 1D理想流体Euler方程xupxeutexpxuutuxuxut1pe11为使激波解稳定，可在压力中加人为粘性项00 0)21(2uuuucqnisnini/pcsniniuuu212116 . 0计算实例：激波管问题20. 高级算法之PIC法F. H. Harlow. 1964, Methods in Comput. Phys., 3, 319发展了PIC (Particle-in-cell) 方法。它综合了Euler格式和Lagrange格式的优点，在Euler网格上定义场变量，用质点表征流体。主要思想是在Euler网格中配置离散的Lagrange流体质点。Euler格式：适

7、用于求解含大扭曲、大滑移问题Lagrange格式：适用于求解撞击界面的多维流体问题。要点：各物理量（除质量外）均定义在网格中心，混合网格 的质量和能量分别按物质给出。初始时质点分布初始时质点分布：以2维均匀网格为例，每个网格内质点数目反映局地密度，且至少要有100个质点；应使计算末每个网格至少要有34个质点（也许应更多）。经验表明，一个网格内质点的分布不宜十分规则，采用菱形分布或不规则分布利于减少质点流过网格边界引起的非物理计算波动。以二维为例，则网格(i, j)的性质为其中的各质点的总和：vvvvvvpeteptt)(1)(0)(KknknknjiKknknjijiyxmjiyxm1,1,)

8、,(),(rvvr及状态方程求出由njinjie,考虑人为粘性项 qvvvvvvqpeteqptt)()(1)(0)(令总共有K个质点，每个质点具有其质量质量 ，坐标，坐标 ，速度速度nkvnkrjyyixxjinknknk)int()int(),(r其中m计算步骤计算步骤1. Euler步以轴对称为例zVqprrUrqpzeVreUteqpzzVVrVUtVqprzUVrUUtU)()(1)()()(动量忽略输运项，上述方程组变为zqVzqVzVprqUrqrUrrrUrpteqpztVqprtU)(1)()()(对动量方程采用FTCS (forward time-centered spac

9、e) 格式求中间速度)(1)(121,21,21,21,21,21,21,21,njinjinjinjijijijinjinjinjinjijijijiqqppztVVqqpprtUU能量网格界面 p 按算术平均值计算，若旁边网格是空网格，置界面 p=0；若界面是固壁，则取界面p等于网格p 。取人为粘性00 0)(22, 1,21uuuqnjinjinjininiUUu1对能量方程采用交错（拉链）差分z向同理)(21)(, 1, 1,21jijijijijiUfUffU若旁边网格是空网格，则U及V均取本网格中心值，若界面是固壁，则法向速度取为0。PIC方法的特点就是利用质点的运动来计算通过网格

10、边界的输运量 。411UaUnk2. 质点步（1）确定质点速度在Euler步中得出的速度 是网格中心的速度，而质点一般不在中心，每个质点的速度可取为与其相邻的四个网格速度的加权平均值（第k个质点）:jijiVU,、设四个网格的公共交点是),(*zrzzzzrrrrnknk*，令)21)(21( )21)(21()21)(21( )21)(21(4321rzarzarzarza4132411VaVnk约定：把覆盖(i, j)网格的部分记为4，对角部分记为1， 沿r方向的部分记为2，沿z方向的部分记为3。4132413241324132（2）计算质点的新坐标tVVzztUUrrnknknknknk

11、nknknk)(21 ,)(211111若覆盖的网格为空网格，空网格的速度取为质点所在网格的速度，或取此网格的面积为0。若覆盖的网格在固壁外，其速度取为相邻内网格的反射速度。（3）重新分配质点质点运动后可能留在原网格，也可能移动到新网格，从而引起网格质量、动量和能量的变化及状态方程求出由1,1,njinjieNknknknjiNknknjijizrmjizrm1111,111,),(),(rvvr缺点对下列情况效果不是非常好：非常不均匀的粒子分布correlated systemscomplex geometries21. 高级算法之SPH法PIC方法同时使用Eulerian和Lagrangi

12、an elements，其欧拉网格的平均密度等于网格中的质点数目。能否不需网格而得到局地密度？ Mesh-freeLucy, 1977, AJ, 82, 1013 Gingold & Monaghan, 1977, MNRAS, 181, 375发展了SPH (Smoothed Particle Hydrodynamics光滑粒子流体动力学) 方法: load each particle with a spatially extended interpolation kernel，平均局地密度缘于所有质点的贡献。自在天体物理研究领域发展以来，已广泛应用于其它物理问题，如连续体结构的解体

13、、碎裂、固体的层裂、脆性断裂等。它可以解决很多其它算法解决不了的问题。由于它是无网格的，故可用于研究极不规则的结构。令 表示为the interpolation kernel centered around 则 r 处的局地密度为（假设总共有N个质点） )(iwrrirNiiiwm1)()(rrrGenerally, a property A(r) is represented by its “smoothed particle estimate”NiiiiiAwmA1)()()()(rrrrrForm of the interpolation kernel: Gaussian or poly

14、nomial. For example:22/32/31)(dsedwsWith a width d chosen such that the number of particles within d is about 5 in 2 dimensions and 15 in 3 dimensions.对动量方程，Monaghan, 1992, ARAA, 30, 543指出直接使用和p会破坏动量和角动量的守恒，并建议使用等式：vvvvvvpeteptt)(1)(0)(21pppNkkkkkkpwmp1)(/rrNiiiwm1)()(rrrikiNkiikkkiwppmt122)(ddv)()(

15、ikikikwwwrrrvvvvvvpeteptt)(1)(0)(ikNkikikkiiwmt1ddvvrikikvvvikiNkikkiwmt1ddv2/ )(kiik其中但是，密度方程不必解此求和，只需由速度计算所有质点的坐标即可：iitvrdd或更准确的表达式其中此关系式满足动量及角动量守恒，且具有相邻质点速度相似的优点。(Monaghan, 1989， JCP, 82, 1)简单的方法是vvvvvvpeteptt)(1)(0)(所有这些方程的求解可采用适当的格式，如蛙跳，龙格库塔等。ikiikNkiikkkiwppmtev122)(dd此外，Gingold & Monaghan, 1977, MNRAS, 181, 375亦考虑了自引力kkkurr kNkukkkuuduududuNGMk1022222/12d)/exp(1)/exp()1(2Example: 中间部变量技巧(Monaghan, 1989， JCP, 82, 1)给定 tn