1.3.1 计算固体力学方法

1.3.1 计算固体力学方法

计算固体力学(computational solid mechanics,CSD)是在已建立的描述固体运动和变形的数学模型的基础上,采用一定的离散化数值方法,用有限个未知量去近似待求的连续函数,从而将微分方程问题转化为代数方程问题,并利用计算机求解[11]。常用的方法包括有限元方法、分子动力学模拟方法、边界元法、无网格方法等。心脑血管植介入器械如血管支架的疲劳性能分析、支架植入过程的仿真和支架与血管壁的相互作用,都是利用计算固体力学的方法进行仿真。

1)有限元方法

有限元方法(finite element method,FEM)是一种在工程中常用的求解偏微分方程的数值方法,其主要思想为将原本连续的结构离散为有限个小单元,将结构上所受到的约束、载荷也离散到对应的节点上,然后利用变分原理进行迭代求解。有限元方法在力学、电磁学、声学、热学等领域的数值仿真上得到了广泛的应用。有限元方法能够求解非常复杂的几何形状和结构的各种力学问题,包括结构静力学问题、结构动力学问题及流体动力学问题等[12]

有限元方法的主要步骤如下:

(1)求解区域的离散化。将连续的结构体剖分为由点、线或面等构成的有限个单元。一般而言,单元的形状可以是任意的,平面问题中常使用三角形、矩形或任意四边形单元等;空间问题中常使用四面体、长方体或任意六面体单元。

(2)选择位移模式。单元的位移模式即位移函数,它是利用单元节点处的位移进行插值得到的,常采用多项式函数作为基函数进行位移的插值。位移模式的选择通常基于单元的类型和单元节点数量。

(3)刚度方程的推导。在确定单元类型和位移模式之后,便可利用最小势能原理等方法建立单元的刚度方程。再将所有单元的刚度方程进行重组,得到整个求解域的刚度方程。

(4)刚度方程的求解。在求解刚度方程时,要充分考虑边界条件的约束,否则可能出现刚度矩阵奇异的结果,进而无法获得数值解。一般而言,刚度矩阵是一个对称、正定、稀疏、带状的矩阵,并且里面含有大量的零元,非零元主要集中于主对角线附近。刚度方程求解的常用方式有直接法和迭代法。直接法中主要是Gauss消去法等,迭代法主要有牛顿-拉夫逊方法(Newton-Raphson method)、高斯-赛德尔迭代法(Gauss-Seidel method)等。

(5)后处理。根据节点位移计算单元的应变和应力,根据实际工况进行分析。(https://www.daowen.com)

有限元方法将连续体进行离散,使用有限个自由度来描述无穷个自由度的系统,相当于提高了原系统的刚度,因而计算所得的位移值总体上将呈现偏小的情况。一般常采用加密网格或使用高阶单元两种方式提高计算精度,前者因为单元构造简单,基函数阶次较低,因而数值稳定性和可靠性更好,但是其收敛性不如后者。虽然如此,在大多数的工程应用中,有限元方法都能够获得符合精度要求的数值解。

2)分子动力学模拟方法

分子动力学模拟(molecular dynamic simulation,MDS)由Alder和Wainwright于1957年首先提出,随后得到快速发展,被广泛应用于材料、能源、化工、生物、医学、激光、电子等众多领域。分子动力学模拟是以分子为基本研究对象,将系统看作具有一定特征的分子集合,运用经典力学或量子力学方法,通过研究微观分子的运动规律,获得各个分子的运动状态,得到体系的宏观特性和基本规律的方法[13-14]。分子动力学模拟是一种借助经典多体体系进行平衡和传递性质计算及过程和现象预测的数值研究方法,对研究对象施加一定的限制条件,随后让其演化到所需要的状态,最后进行测定并提取结果。测定次数越多,结果相应就越准确。

分子动力学模拟以一个包含N个粒子(原子或分子)的体系为对象,在给定粒子之间的作用势、初始条件和边界条件后,通过对牛顿运动方程进行数值积分得到粒子运动轨迹,最后进行统计分析,从微观量获得宏观量。其基本思路如下:

(1)将物质看成由分子和原子组成的粒子系统,从该体系的某一假定的位能模型出发,基于经典力学或量子力学描述的运动规律,再结合粒子的受力情况求解系统中所有粒子在相空间的轨道,进而统计得到系统的热力学参数、结构和输运特性等宏观性质。

(2)根据实际系统所处的不同宏观条件采取不同的系综。常见的系综如下:微正则系综(N,V,E)——孤立、保守系统,与外界无能量交换,体系的粒子数N守恒,体系的体积V保持不变;正则系综(N,V,T)——粒子数N、系统体积V、体系温度T都保持不变,系统总动量为0;等温等压系综(N,P,T)——粒子数N、系统压力P、体系温度T都保持不变。

(3)根据不同的力场(即粒子间势能表达式),通过势能梯度计算粒子间的相互作用力。对于简单粒子系统,常采用Lennard-Jones势、硬球势、软球势、方阱势等模型。

(4)粒子的初始位置通常是随机分布于模拟空间或者由其结晶形态的位置分布决定,初始速度可通过麦克斯韦-玻尔兹曼分布给定。体系常采用周期性、固定、全反射等几类边界条件,具体选择取决于系统所处的环境。

(5)进行趋衡操作以达到目标平衡态,即从系统中增加或移出能量,直到系统能持续给出确定的能量值,此时系统达到平衡。而后即可通过统计的方法计算体系的宏观物理量。