一种土体大变形流动的计算模拟方法
【专利摘要】本发明公开了一种土体大变形流动的计算模拟方法，该方法基于分子动力学（简称MD）方法，以牛顿运动方程为控制方程，引入Hertzian型摩擦公式和斯多克公式分别描述固?固颗粒间相互作用和固?液颗粒间相互作用，运用Verlet算法求解积分过程，建立了固?液耦合的土体大变形流动计算模型。本发明可以呈现土体颗粒复杂流动的基本动力学特征及其运动规律，有利于进一步完善和提高我国防治土体大变形流动灾害的研究水平，同时也为分子动力学计算方法在土体大变形流动灾害中的广泛应用和推广奠定一定基础。
【专利说明】
一种土体大变形流动的计算模拟方法
技术领域
[0001] 本发明属于岩土力学计算领域，涉及一种土体大变形流动的计算模拟方法。
【背景技术】
[0002] 土体大变形流动灾害主要表现为边坡失稳、基坑变形破坏和地面塌陷等诸多形 式，严重威胁人类生命与财产安全。土体大变形的分析对于工程设计和防灾减灾具有重要 的意义，但目前还没有得到较好的解决，其中一个关键原因是数值计算并没有完全反映土 体大变形的力学特征。土体在较大应力状态下的变形过程属于非线性的大变形状态，其变 形不再是符合小变形理论，鉴于传统数值计算方法无法模拟土体大变形的完整过程。因此， 寻求一种有效的计算方法来模拟土体大变形流动过程显得尤为重要。
【发明内容】
[0003] 本发明的目的在于解决现有数值计算方法无法有效模拟土体大变形流动灾害的 技术问题，进而提出一种土体大变形流动的计算模拟方法，本发明将Hertz ian型摩擦公式 和斯多克公式引入分子动力学计算框架，运用Verlet算法求解积分过程，建立了固-液耦合 的土体大变形流动计算模型，可实现土体大变形流动过程的模拟。
[0004] 为达到上述目的，本发明采用如下技术方案：
[0005] -种土体大变形流动的计算模拟方法，其特征在于：所述土体大变形流动的计算 模拟方法包括以下步骤：
[0006] 1)在分子动力学数值计算模拟软件LAMPS中输入粒子信息文件：
[0007] 基于计算体系特征，首先设置计算体系的维数、单位、粒子类型、粒子半径、边界条 件以及积分算法，创建模拟单元;然后设置计算体系基本信息，所述计算体系基本信息包括 但不限于系综、输出的时间间隔、步长以及运行步数;所述计算体系是要进行数值计算模拟 的对象，所述要进行数值计算模拟的对象是土体;所述计算体系特征是要进行数值计算模 拟的对象的物理力学性质及所处边界条件；
[0008] 2)初始设置模块启动，分子动力学数值计算模拟软件LAMMPS对计算体系的所有粒 子赋予初始位置以及初始速度，进行初始化;输入模型参数、粒子初始密度及粒子所受外 力;所述粒子是组成要进行数值计算模拟的对象的分子、原子以及离子；
[0009] 3)模拟控制模块启动，质点邻近搜索：
[0010]根据粒子影响半径，运用Verlet邻近搜索法对计算粒子影响范围内所有粒子进行 搜索，逐一判断周围粒子与计算粒子之间的间距与影响半径之间的关系，确定周围粒子与 计算粒子之间的间距小于影响半径的周围粒子作为计算粒子的邻近粒子;直至确定每个粒 子的所有邻近粒子；
[0011]所述计算粒子是正在被进行计算的粒子；
[0012]所述周围粒子是计算粒子周围的粒子；
[0013]所述影响半径是依据要进行数值计算模拟的对象的颗粒粒径大小设置的距离；
[0014] 4)计算每个粒子的所有邻近粒子对这些粒子所产生的作用力、瞬时加速度以及瞬 时速度：
[0015] 根据步骤2)输入的模型参数、粒子初始密度及粒子所受外力计算步骤3)所得到的 每个粒子的所有邻近粒子对这些粒子所产生的作用力，然后通过矢量求和求得每个粒子上 所受合力，由牛顿运动方程求得瞬时加速度，采用Verlet积分算法对牛顿运动方程求解，得 到每个粒子的瞬时速度和位置；
[0016] 所述每个粒子的所有邻近粒子对这些粒子所产生的作用力的计算方式是采用 Hertzian型摩擦公式和斯多克定律计算得到的，
[0017] 所述Hertzian型摩擦公式是：
[0018]
[0019]其中：
[0020] δ是第i个粒子和第j个粒子之间重叠的距离；
[0021] RhRj分别是第i个粒子以及第j个粒子的半径；
[0022] kn、kt分别是第i个粒子和第j个粒子之间的法向弹性系数以及切向弹性系数;所述 kt = 2/7kn；
[0023] 丫"、丫*分别是第i个粒子和第j个粒子之间的法向阻尼系数以及切向阻尼系数;所 述 y t= 1/2 γ n;
[0024] nieff是第i个粒子和第j个粒子之间的有效质量，所述meff=mimj/(mi+mj);
[0025] mi，mj分别为第i个粒子和第j个粒子的质量；
[0026] Δ St是第i个粒子和第j个粒子之间的切向位移矢量；
[0027] 是连接第i个粒子和第j个粒子中心的单位矢量；
[0028] Vn、Vt分别是第i个粒子和第j个粒子相对速度的法向分量以及切向分量；
[0029] t是时间；
[0030] i是第i个粒子；
[0031] j是除第i个粒子外的所有粒子；
[0032]所述斯多克定律的计算公式是：
[0033] f = -3πη?ν
[0034] 其中：
[0035] η是水分子的粘度；
[0036] d是第i个粒子的粒径；
[0037] V是第i个粒子的运动速度；
[0038] 所述每个粒子的瞬时速度的计算公式是：
[0039]
[0040] 其中：
[0041] i是第i个粒子，即正在被进行计算的粒子；
[0042] t是时间；
[0043] Vi(t)是第i个粒子在t时刻的速度；
[0044] ri(t)是第i个粒子在t时刻的位置；
[0045] At是时间间隔；
[0046] ri(t_ Δ t))是第i个粒子在t_ Δ t时刻的位置；
[0047] ri(t+A t)是第i个粒子在t+Δ t时刻的位置；
[0048] 0(( At)2)是At的二阶无穷小；
[0049] 5)重复步骤4)，对步骤4)所得到的每个粒子的所有邻近粒子对这些粒子所产生的 作用力、瞬时加速度以及瞬时速度进行更新；
[0050] 6)实时计算结果输出模块启动，统计所需物理量，求出模拟系统中粒子在相空间 的轨迹后，应用统计物理原理得到系统的宏观特性和结构特点；
[0051] 7)判断模拟时间是否大于设置的总时间，若模拟时间大于设置的总时间时，步骤 5)的更新过程结束，程序终止;所模拟时间小于设置的总时间时，重复步骤5)直至模拟时间 大于设置的总时间；所述模拟时间是步骤1)中步长与运行步数的乘积;所述设置的总时间 是步骤1)中步长与运行步数总和的乘积。
[0052]本发明的优点如下：
[0053] (1)鉴于目前土体大变形流动问题数值模拟研究的现状，以及传统数值计算方法 的局限性，本发明提供了一种不带近似、动态跟踪粒子轨迹、模拟结果精确，并且具有沟通 宏观特性与微观结构作用的计算模拟方法一一分子动力学方法，能够有效再现土体大变形 流动过程，解决了土体颗粒流动带来的相关轨迹动态变化以及大变形等复杂问题，可为防 治土体大变形流动灾害提供相关指导依据。
[0054] (2)本发明首次将Hertzian型摩擦公式和斯多克公式引入MD框架中，建立了固-液 耦合的土体大变形流动计算模型，该计算模拟方法有效弥补了当前土体大变形流动模拟精 确度不足的问题。
[0055] (3)本发明首次应用分子动力学方法计算模拟土体大变形流动问题。分子动力学 方法不要求模型过分简化，还能克服网格划分带来的问题，提高了计算效率。
【附图说明】
[0056]图1是本发明的流程图；
[0057]图2是砂土流动最终构型的计算模拟结果。
【具体实施方式】
[0058]下面结合附图和具体的实施例对本发明的土体大变形流动的计算模拟方法做进 一步的详细说明：
[0059]参见图1以及图2,本发明提供了一种土体大变形流动的计算模拟方法，该方法包 括以下步骤：
[0060] 1)在分子动力学数值计算模拟软件LAMPS中输入粒子信息文件：
[0061] 基于计算体系特征，首先设置计算体系的维数、单位、粒子类型、粒子半径、边界条 件以及积分算法，创建模拟单元;然后设置计算体系基本信息，计算体系基本信息包括但不 限于系综、输出的时间间隔、步长以及运行步数;计算体系是要进行数值计算模拟的对象， 要进行数值计算模拟的对象是土体;计算体系特征是要进行数值计算模拟的对象的物理力 学性质及所处边界条件；
[0062] 2)初始设置模块启动，分子动力学数值计算模拟软件LAMMPS对计算体系的所有粒 子赋予初始位置以及初始速度，进行初始化;输入模型参数、粒子初始密度及粒子所受外 力;粒子是组成要进行数值计算模拟的对象的分子、原子以及离子；
[0063] 3)模拟控制模块启动，质点邻近搜索：
[0064]根据粒子影响半径，运用Verlet邻近搜索法对计算粒子影响范围内所有粒子进行 搜索，逐一判断周围粒子与计算粒子之间的间距与影响半径之间的关系，确定周围粒子与 计算粒子之间的间距小于影响半径的周围粒子作为计算粒子的邻近粒子;直至确定每个粒 子的所有邻近粒子；
[0065] 计算粒子是正在被进行计算的粒子；
[0066] 周围粒子是计算粒子周围的粒子；
[0067] 影响半径是依据要进行数值计算模拟的对象的颗粒粒径大小设置的距离；
[0068] 4)计算每个粒子的所有邻近粒子对这些粒子所产生的作用力、瞬时加速度以及瞬 时速度：
[0069] 根据步骤2)输入的模型参数、粒子初始密度及粒子所受外力计算步骤3)所得到的 每个粒子的所有邻近粒子对这些粒子所产生的作用力，然后通过矢量求和求得每个粒子上 所受合力，由牛顿运动方程求得瞬时加速度，采用Verlet积分算法对牛顿运动方程求解，得 到每个粒子的瞬时速度和位置；
[0070] 每个粒子的所有邻近粒子对这些粒子所产生的作用力的计算方式是采用 Hertzian型摩擦公式和斯多克定律计算得到的，
[0071 ] Hertzian型摩擦公式是：
[0072]
[0073] 其中：
[0074] δ是第i个粒子和第j个粒子之间重叠的距离；
[0075] RnRj分别是第i个粒子以及第j个粒子的半径；
[0076] kn、kt分别是第i个粒子和第j个粒子之间的法向弹性系数以及切向弹性系数;kt = 2/7kn；
[0077] 丫"、丫*分别是第i个粒子和第j个粒子之间的法向阻尼系数以及切向阻尼系数； T t= 1/2 y η ；
[0078] meff是第i个粒子和第j个粒子之间的有效质量，meff=mimj/(mi+mj);
[0079] nu，!^*别为第i个粒子和第j个粒子的质量；
[0080] △ St是第i个粒子和第j个粒子之间的切向位移矢量；
[0081] 是连接第i个粒子和第j个粒子中心的单位矢量；
[0082] Vn、Vt分别是第i个粒子和第j个粒子相对速度的法向分量以及切向分量；
[0083] t是时间；
[0084] i是第i个粒子；
[0085] j是除第i个粒子外的所有粒子；
[0086]斯多克定律的计算公式：
[0087] f = -3 jrqd v
[0088] 其中：
[0089] η是水分子的粘度；
[0090] d是第i个粒子的粒径；
[0091] V是第i个粒子的运动速度；
[0092] 每个粒子的瞬时速度的计算公式是：
[0093]
[0094] 兵干
[0095] i是第i个粒子，即正在被进行计算的粒子；
[0096] t是时间；
[0097] Vi(t)是第i个粒子在t时刻的速度；
[0098] ri(t)是第i个粒子在t时刻的位置；
[0099] At是时间间隔；
[0100] ri(t_ Δ t))是第i个粒子在t_ Δ t时刻的位置；
[0101 ] r"t- Δ t))是第i个粒子在t+ Δ t时刻的位置；
[0102] 0(( At)2)是At的二阶无穷小；
[0103] 5)重复步骤4)，对步骤4)所得到的每个粒子的所有邻近粒子对这些粒子所产生的 作用力、瞬时加速度以及瞬时速度进行更新；
[0104] 6)实时计算结果输出模块启动，统计所需物理量，求出模拟系统中粒子在相空间 的轨迹后，应用统计物理原理得到系统的宏观特性和结构特点；
[0105] 7)判断模拟时间是否大于设置的总时间，若模拟时间大于设置的总时间时，步骤 5)的更新过程结束，程序终止;所模拟时间小于设置的总时间时，重复步骤5)直至模拟时间 大于设置的总时间；模拟时间是步骤1)中步长与运行步数的乘积;所述设置的总时间是步 骤1)中步长与运行步数总和的乘积。
[0106] 实施例：
[0107]尺寸为40cmX IOcmX 12cm的模型箱水平放置，其中内接箱的有效尺寸为IOcmX IOcmX 10cm，将普通中砂(粒径0.25mm-0.5mm)缓慢均匀地装入内接箱直到8cm高度时停止 装砂，快速抽掉内接箱的挡板使砂土发生流动。
[0108]应用本发明的计算方法，对模型箱内砂土流动进行计算模拟，具体计算步骤如下：
[0109] (1)应用分子动力学软件LAMMPS计算时，设置计算区域长为40cm，宽为10cm，高为 10cm，砂样尺寸长为10cm，宽为10cm，高为8cm，设置体系的单位为约化单位(无量纲），粒子 形状视为圆球形，忽略砂粒的粒径分布，边界条件为固定边界，左右两侧及下侧为wall边 界，上部与外界大气接触，计算直径取粒径期望值的3倍约为2_。
[0110] (2)砂样初始位置与内接箱装好砂土时保持一致，初始速度设为零。
[0111] (3)运用Verlet邻近搜索法，以0.3倍间距，逐一判断周围粒子与计算粒子之间的 间距是否小于影响半径，若小于则该粒子即为邻近粒子，如此往复，逐一确定每个粒子所有 符合条件的邻近粒子，建立邻近列表。
[0112] (4)根据公式（1)和(2)计算体系中粒子间相互作用力，计算参数如表1所示，求出 粒子所受合力，进而求得相应的加速度，再根据公式(3)求得相应的速度。
[0113] (5)重复步骤(4)，逐一更新每一步粒子的位置和速度。
[0114] (6)按照所需输出相应的物理量，如每一步的瞬时速度、加速度、位置、温度等。
[0115] (7)当模拟时间大于设置的总时间时，结束循环计算。最后输出体系的热力学数 据:动能、角动量、体积、密度、时间步和计算速度等。
[0116] 表1计算参数表
[0118]最后说明的是，以上实施例仅用以说明本发明的技术方案而非限制，本领域的普 通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明 技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。
【主权项】
1. 一种土体大变形流动的计算模拟方法，其特征在于:所述土体大变形流动的计算模 拟方法包括以下步骤： 1) 在分子动力学数值计算模拟软件LAMMPS中输入粒子信息文件： 基于计算体系特征，首先设置计算体系的维数、单位、粒子类型、粒子半径、边界条件以 及积分算法，创建模拟单元;然后设置计算体系基本信息，所述计算体系基本信息包括但不 限于系综、输出的时间间隔、步长以及运行步数;所述计算体系是要进行数值计算模拟的对 象，所述要进行数值计算模拟的对象是土体;所述计算体系特征是要进行数值计算模拟的 对象的物理力学性质及所处边界条件； 2) 初始设置模块启动，分子动力学数值计算模拟软件LAMMPS对计算体系的所有粒子赋 予初始位置以及初始速度，进行初始化;输入模型参数、粒子初始密度及粒子所受外力;所 述粒子是组成要进行数值计算模拟的对象的分子、原子以及离子； 3) 模拟控制模块启动，质点邻近搜索： 根据粒子影响半径，运用Verlet邻近搜索法对计算粒子影响范围内所有粒子进行搜 索，逐一判断周围粒子与计算粒子之间的间距与影响半径之间的关系，确定周围粒子与计 算粒子之间的间距小于影响半径的周围粒子作为计算粒子的邻近粒子;直至确定每个粒子 的所有邻近粒子； 所述计算粒子是正在被进行计算的粒子； 所述周围粒子是计算粒子周围的粒子； 所述影响半径是依据要进行数值计算模拟的对象的颗粒粒径大小设置的距离;所述影 响半径SL/2,其中L是模拟单元立方体盒子的边长； 4) 计算每个粒子的所有邻近粒子对这些粒子所产生的作用力、瞬时加速度以及瞬时速 度： 根据步骤2)输入的模型参数、粒子初始密度及粒子所受外力计算步骤3)所得到的每个 粒子的所有邻近粒子对这些粒子所产生的作用力，然后通过矢量求和求得每个粒子上所受 合力，由牛顿运动方程求得瞬时加速度，采用Ver let积分算法对牛顿运动方程求解，得到每 个粒子的瞬时速度和位置； 所述每个粒子的所有邻近粒子对这些粒子所产生的作用力的计算方式是采用 Hertzian型摩擦公式和斯多克定律计算得到的， 所述Hertzian型摩擦公式是：其中： S是第i个粒子和第j个粒子之间重叠的距离； 心、心分别是第i个粒子以及第j个粒子的半径； kn、kt分别是第i个粒子和第j个粒子之间的法向弹性系数以及切向弹性系数;所述kt = 2/7kn； Yn、Yt分别是第i个粒子和第j个粒子之间的法向阻尼系数以及切向阻尼系数;所述 T t= 1/2 y n ； nfeff是第i个粒子和第j个粒子之间的有效质量，所述mrff=mimj/(mi+mj); mi，mj分别为第i个粒子和第j个粒子的质量； A st是第i个粒子和第j个粒子之间的切向位移矢量； 是连接第i个粒子和第j个粒子中心的单位矢量； vn、vt分别是第i个粒子和第j个粒子相对速度的法向分量以及切向分量； t是时间； i是第i个粒子； j是除第i个粒子外的所有粒子； 所述斯多克定律的计算公式是： f = _3Jiqdv 其中： n是水分子的粘度； d是第i个粒子的粒径； v是第i个粒子的运动速度； 所述每个粒子的瞬时速度的计算公式是：其中： i是第i个粒子，即正在被进行计算的粒子； t是时间； Vi(t)是第i个粒子在t时刻的速度； n (t)是第i个粒子在t时刻的位置； At是时间间隔； ri(t_ A t))是第i个粒子在t-A t时刻的位置； ri(t+A t)是第i个粒子在t+A t时刻的位置； 〇(( At)2)是At的二阶无穷小； 5) 重复步骤4)，对步骤4)所得到的每个粒子的所有邻近粒子对这些粒子所产生的作用 力、瞬时加速度以及瞬时速度进行更新； 6) 实时计算结果输出模块启动，统计所需物理量，求出模拟系统中所有粒子在相空间 的轨迹后，应用统计物理原理得到计算体系的宏观特性和结构特点； 7) 判断模拟时间是否大于设置的总时间，若模拟时间大于设置的总时间时，步骤5)的 更新过程结束，程序终止;所模拟时间小于设置的总时间时，重复步骤5)直至模拟时间大于 设置的总时间；所述模拟时间是步骤1)中步长与运行步数的乘积;所述设置的总时间是步 骤1)中步长与运行步数总和的乘积。
【文档编号】G06F17/50GK106055851SQ201610573468
【公开日】2016年10月26日
【申请日】2016年7月19日
【发明人】赵紫阳, 申俊敏, 张军, 孙志杰, 薛晓辉, 马林
【申请人】山西省交通科学研究院, 山西交科公路勘察设计院