非线性方程全隐式数值求解方法及其在地球流体力学中的应用.docxVIP

非线性方程全隐式数值求解方法及其在地球流体力学中的应用 1 全隐式差分格式 灾难报告一直是气象科学的一个重要课题。国际上正着手开展的 “THORPEX (The Observing system Research and Predictability Experiment) ” 计划, 其目的就在于改进对有重大影响 (high-impact) 天气事件的预报。作为主要预报手段, 发展完善数值预报模式就成为当务之急。它涉及多方面的问题, 其中提高模式的水平和垂直分辨率是基本的, 也是有效的手段, 这一点已为近二三十年来数值预报实践所肯定。近年来, 计算机和计算技术迅猛发展, 充分利用计算机资源, 极大地提高模式分辨率, 发展高精度新一代数值模式, 已成为各国数值预报研究的共同努力方向, 不仅是针对中小尺度系统区域模式, 对于中期预报和短期气候预报模式也是如此。例如, 欧洲中期天气预报中心的计划, 拟在2008年全球模式达到15 km分辨率, 而实际进展, 很可能提前实现。 然而提高分辨率, 不是简单地缩小网格就可达到预期目的。研究表明, 在目前模式水平上继续提高分辨率是否有明显收益, 采用全隐式方案是保证之一。因此, 发展全隐式模式, 研究它的高效算法, 成为模式发展的一个前沿课题。 多年来, 各国的模式广泛采用半隐式差分格式, 对不同时间尺度的大气运动采取不同的时间差分格式。从计算稳定性看, 半隐式格式容许人们取时间步长远比显式格式的长。然而, 这种在方程组中对产生快波和慢波项分别采用显式和隐式差分的做法, 造成不协调形成新的误差源。例如, 对σ坐标及本文使用的高度地形追随坐标?zz?中气压梯度项的差分, 其第一项是线性项直接做隐式差分; 而第二项是非线性项, 先线性化得出一项线性项和一项为偏差变量的非线性项, 前者使用隐式差分, 后者使用显式差分。如此在气压梯度项这三部分中, 使用两种不同差分格式, 无疑将产生不协调, 引起计算误差。另外, 很难做到差分格式保持原方程组的时空守恒性。很多采用半隐式格式的天气预报和气候模式中, 所使用的总能量守恒计算格式, 仅是在显式格式条件下才成立。例如欧洲中期天气预报中心 (ECMWF) 原欧拉格式的模式就是如此。只有向全隐式格式发展, 才有助于解决上述问题。 在这组非线性方程中取全隐式格式, 即不管是在线性项还是非线性项中的变量均取下一时刻 (一阶精度) 或取下一和当前时刻的平均 (二阶精度), 就能克服上面的问题。用迭代求解, 在给定的时间步长内多次迭代直至非线性线误差很小, 就称之为非线性协调求解方法。在这种全隐式非线性协调求解方法中就不会出现像半隐式格式那样产生不协调引起的误差及对计算格式守恒性影响 (如果此格式是守恒的)。我国学者早已认识采用全隐式格式的必要, 并沿此方向长期不懈努力, 例如, 有名的全隐式能量守恒格式, 在离散化的差分方程组中应用此格式可使整个系统在时间和空间上保持总能量守恒, 和防止能量向短波单向传递的 “尺度守恒性”, 保证计算绝对稳定。这比世界上常用的瞬时空间守恒格式进了一大步。然而, 这种全隐式格式尚未解决迭代收敛问题, 致使在实际应用时, 不管收敛与否只要达到事先设定迭代次数, 就停止迭代。这会带来不小的误差, 并且总能量难于保持原设计的时空守恒, 因而需要探索更好的求解方法。 此外, 大气谱模式, 是对产生重力波的项进行线性化, 然后作隐式处理, 其他项如平流项作显式处理, 最后消元形成对一个变量的线性Helmholtz或Poisson方程。全隐式格式应用于大气非线性方程组时, 就难于形成一个线性的Helmholtz方程。另外, 对于非静力的半隐式半拉格朗日预报模式, 最后要求解一个三维Helmholtz方程。这个方程是变系数, 且带交叉微商项, 很难用谱方法处理。至今还没有看到用全隐式格式求解谱模式的报导, 这可能是一个原因。 近年来, 在计算数学领域内发展出一种专门针对非线性方程隐式求解方法, 称为JFNK(Jacobin-free Newton-Krylov) 方法。Newton方法提供非线性协调性, 而Krylov方法是用来解大型矩阵。JFNK方法一个突出特点是, 它不要求形成和存贮Jacobian矩阵, 只是将Jacobian矩阵的作用通过一种差分来近似求得。其次, 在JFNK方法中用Krylov子空间方法求解矩阵时往往使用广义的最小余差法 (GMRES), 它是求解大型非对称矩阵的有效方法。十分关键的是,JFNK引入预条件处理器 (preconditioner), 或称 “预条件子”, 使求解的矩阵尽量对角化, 迭代求解十分有效。预条件处理器的有效构造, 十分依赖于所研究的物理问题, 因而具有充分的研究发展空间, 需要我们气象工作者去研