有限元分析快速直接求解技术进展

现代的有限元分析往往产生大规模的线性方程组, 它的求解效率是有限元分析 中最关键的一环. 自20世纪90年代中期, 有限元的求解技术发生了巨大的变化, 传统的变带 宽解法与波前法被稀疏解法所替代. 这一替代为有限元分析带来了求解速度的突破, 它使得1万到10万个节点的实用三维有限 元分析在微机上即时求解成为现实. 本文回顾非并行有限元快速直接求解技术在过去20年 的发展, 着重讨论了填充元优化与浮点加速运算方法, 期望能引起同行的注意.     

基于局部通信的有限元分析并行算法优化研究

针对有限元分析的计算问题,在现有采用全局通信方案的简单并行算法基础上,对其所涉核心算法,采用稀疏数据结构与局部通信进行并行算法优化设计,有效减少了通信所涉及的处理器个数与通信量.同时,通过采用非阻塞通信,并将与通信无关计算进行分离与前置的方法,进行计算与通信重叠,以有效隐藏通信开销的影响.实验结果表明,优化所得算法相比现有算法具有明显改进,特别是对稀疏矩阵稠密向量乘与单元贡献装配,改进很大.同时,随着任务个数的增加,改进效果越明显.     

轴弯曲与不平衡柔性转子共振稳态响应随机分析

系统参数随机性是影响准确预测转子动力学行为的重要因素。为此,研究了具有初始轴弯曲的不平衡柔性转子在故障参数随机情况下的不确定性分析问题。首先,基于转子动力学梁元有限元理论,推导了不平衡和轴弯曲故障共同作用下的系统稳态动力学方程,并以轴心轨迹长半轴作为关键响应量建立了柔性转子共振稳态响应与输入参数间的模型函数;其次,联合广义多项式混沌展开、留一法交叉验证以及最小角回归技术实现了柔性转子共振稳态响应的非嵌入式自适应稀疏多项式混沌展开,并与基于普通最小二乘法的多项式混沌展开和基于Monte Carlo仿真的结果作了对比分析,验证了自适应稀疏展开方案的有效性、精度和效率;最后,以构建的自适应稀疏多项式混沌展开式作为近似模型,重点分析了转子圆盘处一阶共振稳态响应的随机特性,并基于Sobol指标获得了响应对各故障参数的全局灵敏度指标。     

混凝土细观力学分析程序中的快速算法与并行算法设计

针对一套混凝土细观力学分析程序,在分析其计算方法与计算效率的不足之后,提出了采用稀疏矩阵与稀疏向量技术来高效实现有限元刚度矩阵装配过程的算法,并采用双门槛不完全Cholesky分解预条件技术与CG法相结合来高效地求解稀疏线性方程组。之后,从整体上提出了一个将有限单元分布与未知量分布有机结合的并行算法设计方案,并分别针对刚度矩阵装配、双门槛不完全Cholesky分解、稀疏矩阵与稠密向量相乘、稀疏向量相加等核心算法,进行了相应的并行算法设计。最后,在由每节点2 CPU的8个Intel Xeon节点采用千兆以太网连成的机群上,针对两个混凝土数值试样进行了数值实验,第一个试样含44117个网格点与53200个有限单元,第二个试样含71013个网格点与78800个有限单元;对第一个试样,原串行程序进行全程567次加载计算需要984.83小时约41天,采用文中串行算法后,模拟时间减少到22531秒约6.26小时,采用并行算法在16个CPU上的模拟时间进一步降为3860秒约1.07小时。对第二个试样,原串行程序进行全程94次加载计算需要467.19小时约19.5天,采用文中串行算法后,模拟时间减少到11453秒约3.18小时,采用并行算法在16个CPU上的模拟时间进一步降为1704秒约28.4分钟。串行算法的改进与并行算法的设计大大缩短了计算时间,对加快混凝土力学性能的分析研究具有重要意义。     

金属结构航天器陨落过程三维瞬态传热有限元算法研究

构建金属桁架结构航天器陨落再入气动热环境有限元传热模型，是准确预测在轨服役期满大型航天器陨落再入解体过程温度分布的关键。本文采用四节点四面体单元对空间进行离散，依据泛函理论，将传热控制方程离散为代数方程组；利用有限单元法总体合成得到具有对称正定、高度稀疏和非0元素分布的规则性刚度矩阵，发展一维变带宽压缩存贮技术，有效解决大型稀疏矩阵的数据存贮问题；为有效抑制求解过程出现的温度在时间和空间上的振荡问题，发展集中热容矩阵系数处理方法，将热容矩阵的同行或同列元素相加代替对角线元素，使非对角线元素化为0，构造求解三维瞬态温度场的两点向后差分格式、Crank-Nicolson格式和Galerkin格式。通过对正方体瞬态传热计算验证分析，在相同条件下，采用以上三种格式均可获得一致稳定的温度解，并得到与现有ANSYS有限元软件较为吻合的计算结果，验证了所建立三维瞬态传热有限元计算模型的准确可靠性。在此基础上，对铝合金低轨航天器薄壳结构进行了传热计算，给出了类天宫飞行器两舱体陨落飞行107.5 km~90 km不同高度的瞬态温度分布，为这类寿命末期航天器陨落再入解体预报提供理论支撑与可计算模型。     

侧向稀疏对PZT—95／5陶瓷冲击去极化影响的研究

采用平面冲击波加载，对样品在不同封装介绍（不同侧向稀疏程度的情况下）和几何尺寸条件下进行冲击去极化实验，结果表明侧向稀疏导致电荷释放以近乎线性关系减小；根据实验条件拟合出物理计算模型，利用流体弹塑性模型，采用有限元计算方法对样品中的冲击压力进行三维计算，计算结果能很好地解释实验现象。     

一种新的运动激波--非定常边界层相互干扰现象

本文探讨了一种新的激波－非定常边界层相互干扰现象，这种激波－边界层干扰现象既不同于定常激波－边界层干扰现象，又不同于激波在端面反射后与该激波所诱导的边界层之间的干扰现象，而是运动激波与稀疏波和第一激波所诱导的非这常边界层之间的干扰现象，本文对这种现象用微波动力学理论进行分析，并把这种干扰现象看成激波的绕射现象，同时在稀疏波破膜的双驱动激波管中进行实验观察，最后把理论分析与实验观察进行了比较。     

精密离心机结构动静态特性计算

精密离心机是测量和标定加速度计的重要设备，而加速度计的精度直接关系到运载体的导航精度，本文结合精密离心机的机械结构设计，利用结构分析的有限元方法，对离心机结构设计的方案进行了详细的结构动静态特性分析计算，得出了精密离心机的动静态特性。     

一种有限元模型动力缩聚移频迭代法

提出了一种基于矩阵广义逆的有限元模型动力缩聚移频迭代方法,该方法首先直接从原系统特征方程出发,导出反映系统主,副自由度之间位移关系的动力缩聚矩阵的控制方程,然后给出了相应的迭代求解方法和收敛准则。为了减少求矩阵广义逆的计算工作量,本文给出了一种替代方法,把对一个高阶满阵求逆转化为对一个同阶高度稀疏矩阵求逆。与已有的动力缩聚迭代法相比,本文提出的方法具有两个显著的优点：其一是迭代收敛速度高,其二是通     

一种求解有限元问题结点平衡方程的快速方法

目前随着工程实际问题复杂程度的增加及分析的要求,特别是材料非线性分析的引入,尽管计算机的运算速度、内存、外存容量等不断提高,但并不能完全满足大规模计算的需要,更快、更节省存贮空间的算法一直是有限元法分析过程中的一项核心技术要求。本文针对一些结点及单元均规律化地排列的有限元问题,提出可以透过其相邻结点的关系记录结点平衡方程中系数矩阵的非零元素,无需再像等带宽存贮那样去记录带宽内大量的零元素。此方法可以大大地减少系数矩阵元素的存贮量,从而可以提高计算机读取数据的速度及改善利用迭代法求解的效率。     

A NEW HIGH PERFORMANCE SPARSE STATIC SOLVER IN FINITE ELEMENT ANALYSIS WITH LOOP-UNROLLING

In the previous papers, a high performance sparse static solver with two-level unrolling based on a cell-sparse storage scheme was reported. Although the solver reaches quite a high efficiency for a big percentage of finite element analysis benchmark tests, the MFLOPS （million floating operations per second） of LDL^T factorization of benchmark tests vary on a Dell Pentium IV 850 MHz machine from 100 to 456 depending on the average size of the super-equations, i.e., on the average depth of unrolling. In this paper, a new sparse static solver with two-level unrolling that employs the concept of master-equations and searches for an appropriate depths of unrolling is proposed. The new solver provides higher MFLOPS for LDL^T factorization of benchmark tests, and therefore speeds up the solution process.     

Efficient approximation of Sparse Jacobians for time‐implicit reduced order models

This paper introduces a sparse matrix discrete interpolation method to effectively compute matrix approximations in the reduced order modeling framework. The sparse algorithm developed herein relies on the discrete empirical interpolation method and uses only samples of the nonzero entries of the matrix series. The proposed approach can approximate very large matrices, unlike the current matrix discrete empirical interpolation method, which is limited by its large computational memory requirements. The empirical interpolation indices obtained by the sparse algorithm slightly differ from the ones computed by the matrix discrete empirical interpolation method as a consequence of the singular vectors round‐off errors introduced by the economy or full singular value decomposition (SVD) algorithms when applied to the full matrix snapshots. When appropriately padded with zeros, the economy SVD factorization of the nonzero elements of the snapshots matrix is a valid economy SVD for the full snapshots matrix. Numerical experiments are performed with the 1D Burgers and 2D shallow water equations test problems where the quadratic reduced nonlinearities are computed via tensorial calculus. The sparse matrix approximation strategy is compared against five existing methods for computing reduced Jacobians: (i) matrix discrete empirical interpolation method, (ii) discrete empirical interpolation method, (iii) tensorial calculus, (iv) full Jacobian projection onto the reduced basis subspace, and (v) directional derivatives of the model along the reduced basis functions. The sparse matrix method outperforms all other algorithms. The use of traditional matrix discrete empirical interpolation method is not possible for very large dimensions because of its excessive memory requirements. Copyright © 2016 John Wiley & Sons, Ltd.     

Application of conjugate-gradient-like methods to a hyperbolic problem in porous-media flow

This paper presents the application of a preconditioned conjugate-gradient-like method to a non-self-adjoint problem of interest in underground flow simulation. The method furnishes a reliable iterative solution scheme for the non-symmetric matrices arising at each iteration of the non-linear time-stepping scheme. The method employs a generalized conjugate residual scheme with nested factorization as a preconditioner. Model runs demonstrate significant computational savings over direct sparse matrix solvers.     

一种简单的局部离散速度统一气体动理论格式

统一气体动理论格式UGKS(Unified Gas-Kinetic Scheme)是一种适用于从连续流到自由分子流的全流域计算格式。在该格式中一般使用统一的离散速度空间。而在高速流动中,不同节点的分布函数往往差异很大。为了保证计算的精度,离散速度空间必须满足所有节点的需要,占用了大量的内存。采用局部的均匀离散速度空间,离散速度的范围随节点状态的变化而变化,从而降低了内存的需要,并通过引入背景网格避免了不同节点离散速度的插值。最后,通过两个一维算例对该方法进行了测试。测试结果显示,采用局部离散速度空间能够得到可靠的结果,并且在模拟高速流动时计算效率明显提高。     

A hybrid Padé ADI scheme of higher‐order for convection–diffusion problems

A high‐order Padé alternating direction implicit (ADI) scheme is proposed for solving unsteady convection–diffusion problems. The scheme employs standard high‐order Padé approximations for spatial first and second derivatives in the convection‐diffusion equation. Linear multistep (LM) methods combined with the approximate factorization introduced by Beam and Warming (J. Comput. Phys. 1976; 22 : 87–110) are applied for the time integration. The approximate factorization imposes a second‐order temporal accuracy limitation on the ADI scheme independent of the accuracy of the LM method chosen for the time integration. To achieve a higher‐order temporal accuracy, we introduce a correction term that reduces the splitting error. The resulting scheme is carried out by repeatedly solving a series of pentadiagonal linear systems producing a computationally cost effective solver. The effects of the approximate factorization and the correction term on the stability of the scheme are examined. A modified wave number analysis is performed to examine the dispersive and dissipative properties of the scheme. In contrast to the HOC‐based schemes in which the phase and amplitude characteristics of a solution are altered by the variation of cell Reynolds number, the present scheme retains the characteristics of the modified wave numbers for spatial derivatives regardless of the magnitude of cell Reynolds number. The superiority of the proposed scheme compared with other high‐order ADI schemes for solving unsteady convection‐diffusion problems is discussed. A comparison of different time discretizations based on LM methods is given. Copyright © 2009 John Wiley & Sons, Ltd.     

Effects of the Jacobian evaluation on Newton's solution of the Euler equations

Newton's method is developed for solving the 2‐D Euler equations. The Euler equations are discretized using a finite‐volume method with upwind flux splitting schemes. Both analytical and numerical methods are used for Jacobian calculations. Although the numerical method has the advantage of keeping the Jacobian consistent with the numerical residual vector and avoiding extremely complex analytical differentiations, it may have accuracy problems and need longer execution time. In order to improve the accuracy of numerical Jacobians, detailed error analyses are performed. Results show that the finite‐difference perturbation magnitude and computer precision are the most important parameters that affect the accuracy of numerical Jacobians. A method is developed for calculating an optimal perturbation magnitude that can minimize the error in numerical Jacobians. The accuracy of the numerical Jacobians is improved significantly by using the optimal perturbation magnitude. The effects of the accuracy of numerical Jacobians on the convergence of the flow solver are also investigated. In order to reduce the execution time for numerical Jacobian evaluation, flux vectors with perturbed flow variables are calculated only for neighbouring cells. A sparse matrix solver that is based on LU factorization is used. Effects of different flux splitting methods and higher‐order discretizations on the performance of the solver are analysed. Copyright © 2005 John Wiley & Sons, Ltd.     

A second-order upwind finite-volume method for the Euler solution on unstructured triangular meshes

A scheme for the numerical solution of the two-dimensional (2D) Euler equations on unstructured triangular meshes has been developed. The basic first-order scheme is a cell-centred upwind finite-volume scheme utilizing Roe's approximate Riemann solver. To obtain second-order accuracy, a new gradient based on the weighted average of Barth and Jespersen's three-point support gradient model is used to reconstruct the cell interface values. Characteristic variables in the direction of local pressure gradient are used in the limiter to minimize the numerical oscillation around solution discontinuities. An Approximate LU (ALU) factorization scheme originally developed for structured grid methods is adopted for implicit time integration and shows good convergence characterisitics in the test. To eliminate the data dependency which prohibits vectorization in the inversion process, a black-gray-white colouring and numbering technique on unstructured triangular meshes is developed for the ALU factorization scheme. This results in a high degree of vectorization of the final code. Numerical experiments on transonic Ringleb flow, transonic channel flow with circular bump, supersonic shock reflection flow and subsonic flow over multielement aerofoils are calculated to validate the methodology.     

A high‐order alternating direction implicit method for the unsteady convection‐dominated diffusion problem

A high‐order alternating direction implicit (ADI) method for solving the unsteady convection‐dominated diffusion equation is developed. The fourth‐order Padé scheme is used for the discretization of the convection terms, while the second‐order Padé scheme is used for the diffusion terms. The Crank–Nicolson scheme and ADI factorization are applied for time integration. After ADI factorization, the two‐dimensional problem becomes a sequence of one‐dimensional problems. The solution procedure consists of multiple use of a one‐dimensional tridiagonal matrix algorithm that produces a computationally cost‐effective solver. Von Neumann stability analysis is performed to show that the method is unconditionally stable. An unsteady two‐dimensional problem concerning convection‐dominated propagation of a Gaussian pulse is studied to test its numerical accuracy and compare it to other high‐order ADI methods. The results show that the overall numerical accuracy can reach third or fourth order for the convection‐dominated diffusion equation depending on the magnitude of diffusivity, while the computational cost is much lower than other high‐order numerical methods. Copyright © 2011 John Wiley & Sons, Ltd.     

Running unstructured grid‐based CFD solvers on modern graphics hardware

Techniques used to implement an unstructured grid solver on modern graphics hardware are described. The three‐dimensional Euler equations for inviscid, compressible flow are considered. Effective memory bandwidth is improved by reducing total global memory access and overlapping redundant computation, as well as using an appropriate numbering scheme and data layout. The applicability of per‐block shared memory is also considered. The performance of the solver is demonstrated on two benchmark cases: a NACA0012 wing and a missile. For a variety of mesh sizes, an average speed‐up factor of roughly 9.5 × is observed over the equivalent parallelized OpenMP code running on a quad‐core CPU, and roughly 33 × over the equivalent code running in serial. Copyright © 2010 John Wiley & Sons, Ltd.