CN117610393B - 刚性弹体侵彻多层土体深度的数值仿真预测方法及系统 - Google Patents

刚性弹体侵彻多层土体深度的数值仿真预测方法及系统 Download PDF

Info

Publication number
CN117610393B
CN117610393B CN202311609709.7A CN202311609709A CN117610393B CN 117610393 B CN117610393 B CN 117610393B CN 202311609709 A CN202311609709 A CN 202311609709A CN 117610393 B CN117610393 B CN 117610393B
Authority
CN
China
Prior art keywords
rigid
soil
particles
penetration
projectile
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202311609709.7A
Other languages
English (en)
Other versions
CN117610393A (zh
Inventor
陈丁
黄文雄
黄丹
高景芝
袁德荣
姚学昊
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hohai University HHU
Original Assignee
Hohai University HHU
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Hohai University HHU filed Critical Hohai University HHU
Priority to CN202311609709.7A priority Critical patent/CN117610393B/zh
Publication of CN117610393A publication Critical patent/CN117610393A/zh
Application granted granted Critical
Publication of CN117610393B publication Critical patent/CN117610393B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/25—Design optimisation, verification or simulation using particle-based methods
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
    • G16C60/00—Computational materials science, i.e. ICT specially adapted for investigating the physical or chemical properties of materials or phenomena associated with their design, synthesis, processing, characterisation or utilisation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/26—Composites
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Computing Systems (AREA)
  • General Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Geometry (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Investigating Strength Of Materials By Application Of Mechanical Stress (AREA)

Abstract

本发明公开了刚性弹体侵彻多层土体深度的数值仿真预测方法及系统，属于冲击动力学技术领域。方法包括：S1、构建刚性弹体三维几何模型以及土体分布几何模型，并对刚性弹体三维几何模型以及土体分布几何模型进行参数赋值，以及网格划分；S2、以网格中心为粒子坐标，网格体积为粒子体积，生成刚性弹体粒子离散模型以及土体粒子离散模型；S3、初始化刚性弹体及土体的信息，分别计算土体粒子的运动方程以及刚性弹体粒子的运动方程；S4、计算刚性弹体与土体之间的接触摩擦力，基于接触摩擦力、土体粒子的运动方程以及刚性弹体粒子的运动方程计算刚性弹体的运动轨迹，记录刚性弹体速度为0时的侵彻水平深度和竖直深度。

Description

刚性弹体侵彻多层土体深度的数值仿真预测方法及系统
技术领域
本发明属于冲击动力学技术领域，具体涉及刚性弹体侵彻多层土体深度的数值仿真预测方法及系统。
背景技术
钻地弹侵彻地层的侵彻深度是用来评估地下工程和国防安全的重要参数。当侵彻深度不同时，地下结构物会表现出不同程度的损伤效果。工程上，大多数地下防护工程具有多层材料、深埋等防护特点。其中，地下结构物上层介质主要是以土体为主。土体又是由颗粒材料组成，而颗粒材料既有固体的承载力特性，也有流体的流动特性。由于土体的多样性、离散性，实际地层大多数呈现层状，且每一层的土体材料属性均不相同。同时，弹体本身的三维运动姿态也是主要影响侵彻深度的关键因素。通过精确评估弹体的多层土体侵彻深度，可对地下工程结构的安全防护设计提供设计指导。
侵彻深度受到土体材料属性、弹体的不同形态及姿势、弹体与土体材料之间的界面接触模型等重要因素影响。弹体侵彻多层土体的运动规律主要是通过试验、理论分析以及数值仿真得到。试验主要是通过对试验数据进行拟合，但试验价格较昂贵，难以大量重复试验，且不能很好地涵盖复杂多层土体的多样性。理论分析存在简化过多、对于弹体三维运动状态的评估不准确等问题，大多数被用来预测垂直分层土体的弹体正向侵彻。数值仿真可以将弹体简化为刚体，并设定三维运动形态，且可以很好地描述复杂的多层土体以及其土体强度对弹体侵彻深度的影响。但传统网格类方法在数值仿真计算过程中土体网格模型容易出现网格畸变，其预测精度和有效性难以保证。同时，现有数值仿真技术大多数并未考虑弹体和土体的界面接触摩擦。因此，上述手段均较难对刚体弹体侵彻多层土体的深度进行预测，不能综合考虑多种关键参数对弹体侵彻深度的影响，并做出正确评价。
发明内容
本发明旨在解决现有技术的不足，提出一种刚性弹体侵彻多层土体的侵彻深度数值仿真预测方法及系统，考虑到光滑粒子流体动力学方法在大变形问题中的优势，对刚性弹体侵彻土体过程进行光滑粒子流体动力学方法(SPH)三维模拟，通过耦合土体强度以及刚性弹体和土体摩擦接触界面作用这两个主要的影响侵彻深度因素，突破了传统方法网格畸变问题，综合考虑了不同形态及姿势的刚性弹体、土体材料属性、刚性弹体与土体材料之间的界面接触摩擦作用。
为实现上述目的，本发明提供了如下方案：一种刚性弹体侵彻多层土体深度的数值仿真预测方法，包括以下步骤：
S1、构建刚性弹体三维几何模型以及土体分布几何模型，并对所述刚性弹体三维几何模型以及土体分布几何模型进行参数赋值，以及网格划分；
S2、以网格中心为粒子坐标，网格体积为粒子体积，生成刚性弹体粒子离散模型以及土体粒子离散模型；
S3、初始化刚性弹体及土体的信息，分别计算土体粒子的运动方程以及刚性弹体粒子的运动方程；
S4、计算刚性弹体与土体之间的接触摩擦力，基于所述接触摩擦力、所述土体粒子的运动方程以及所述刚性弹体粒子的运动方程计算刚性弹体的运动轨迹，记录刚性弹体速度为0时的侵彻水平深度和竖直深度。
进一步优选地，所述土体粒子的运动方程包括：
式中，表示为时间微分；ma、ρa分别表示土体粒子a的质量和质量密度；va表示土体粒子a的速度矢量；/>表示侵彻外力；/>表示内部约束力；gab土体粒子a对土体粒子b表示辅助核函数导数的离散形式；b表示土体粒子b；Ωa表示土体粒子a的支持域；Vb表示土体粒子b的体积；vb表示土体粒子b的速度矢量。
进一步优选地，所述刚性弹体粒子的运动方程包括：
式中，表示为时间微分；M0表示刚性弹体粒子质量；v0表示刚性弹体的速度；Text表示刚性弹体受到的抗侵彻力；G表示刚性弹体粒子重力；J0表示刚性弹体三个方向上的惯性矩；Sext表示刚性弹体受到的抗侵彻外合力矩；ω0表示刚性弹体的三个方向上的角速度。
进一步优选地，所述所述接触摩擦力包括：法向接触力fn和切向接触摩擦力fτ，并分别作用累加到土体粒子受到的侵彻外力、以及刚性弹体受到的抗侵彻力及抗侵彻外合力矩；
所述土体粒子受到的侵彻外力包括：
式中，fn,a表示法向接触力；fτ,a表示切向接触摩擦力；
所述刚性弹体受到的抗侵彻力包括：
式中，a表示土体粒子；Bs表示土体粒子集合；
所述刚性单体受到的抗侵彻外合力矩包括：
式中，xa0表示刚性弹体重心到土体粒子a的距离。
本发明还提供一种刚性弹体侵彻多层土体深度的数值仿真预测系统，所述系统用于实现上述方法，包括：划分单元、生成单元、计算单元及轨迹单元；
所述划分单元用于构建刚性弹体三维几何模型以及土体分布几何模型，并对所述刚性弹体三维几何模型以及土体分布几何模型进行参数赋值，以及网格划分；
所述生成单元用于以网格中心为粒子坐标，网格体积为粒子体积，生成刚性弹体粒子离散模型以及土体粒子离散模型；
所述计算单元用于初始化刚性弹体及土体的信息，分别计算土体粒子的运动方程以及刚性弹体粒子的运动方程；
所述轨迹单元用于计算刚性弹体与土体之间的接触摩擦力，基于所述接触摩擦力、所述土体粒子的运动方程以及所述刚性弹体粒子的运动方程计算刚性弹体的运动轨迹，记录刚性弹体速度为0时的侵彻水平深度和竖直深度。
进一步优选地，所述土体粒子的运动方程包括：
式中，表示为时间微分；ma、ρa分别表示土体粒子a的质量和质量密度；va表示土体粒子a的速度矢量；/>表示侵彻外力；/>表示内部约束力；gab土体粒子a对土体粒子b表示辅助核函数导数的离散形式；b表示土体粒子；Ωa表示土体粒子a的支持域；Vb表示土体粒子b的体积；vb表示土体粒子b的速度矢量。
进一步优选地，所述刚性弹体粒子的运动方程包括：
式中，表示为时间微分；M0表示刚性弹体质量；v0表示刚性弹体的速度；Text表示刚性弹体受到的抗侵彻力；G表示刚性弹体重力；J0表示刚性弹体三个方向上的惯性矩；Sext表示刚性弹体受到的抗侵彻外合力矩；ω0表示刚性弹体的三个方向上的角速度。
进一步优选地，所述接触摩擦力包括：法向接触力fn和切向接触摩擦力fτ，并分别作用累加到土体粒子受到的侵彻外力、以及刚性弹体受到的抗侵彻力及抗侵彻外合力矩；
所述土体粒子受到的侵彻外力包括：
式中，fn,a表示法向接触力；fτ,a表示切向接触摩擦力；
所述刚性弹体受到的抗侵彻力包括：
式中，a表示土体粒子；Bs表示土体粒子集合；
所述刚性单体受到的抗侵彻外合力矩包括：
式中，xa0表示刚性弹体重心到土体粒子a的距离。
与现有技术相比，本发明的有益效果为：
本发明建立了SPH法实现了刚性弹体侵彻多层土体的数值仿真技术，不仅包含了SPH法用于求解土体大变形的优势，同时还引入刚性弹体和土体接触摩擦的点对体积接触摩擦模型，点对体积接触摩擦算法通过土体粒子支持域覆盖刚性弹体体积得到，可以计算刚性弹体侵彻土体过程中界面接触材料参数对于深度的预测影响，从而实现侵彻全过程的模拟。为刚性弹体侵彻复杂多层土体的侵彻深度预测提供新的途径。
本发明不仅克服了传统试验方法中大量重复实施所带来开销大的缺陷，而且避免了理论中过多简化、难以精确预测的不足。建立的刚体侵彻多层土体的过程数值仿真方法可以有效地获取刚性弹体侵彻深度，还可以通过调整不同的三维刚性弹体形态以及材料参数，在低成本下进行大量的模拟仿真。
附图说明
为了更清楚地说明本发明的技术方案，下面对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例刚性弹体侵彻多层土体深度的数值仿真预测方法的流程示意图；
图2为本发明实施例粒子离散以及土体粒子对刚性弹体之间的接触关系示意图；
图3为本发明实施例测试算例的刚性弹体、土体几何尺寸示意图；
图4为本发明实施例初始参数下刚性弹体侵彻的最终侵彻深度结果示意图；
图5为本发明实施例调整参数下刚性弹体侵彻的最终侵彻深度结果示意图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
实施例一：
如图1所示，本实施例提供一种刚性弹体侵彻多层土体深度的数值仿真预测方法，包括以下步骤：
S1、构建刚性弹体三维几何模型以及土体分布几何模型，并对所述刚性弹体三维几何模型以及土体分布几何模型进行参数赋值，以及网格划分。
根据实际情况创建不同分层形状的土体分布几何模型以及刚性弹体三维几何模型；对每一层土体分别赋值不同密度ρ，以及其材料本构参数：弹性模量E、泊松比ν、内摩擦角黏聚力cm以及剪胀角ψ等。设置刚性弹体与土体之间的法向接触刚度kn、切向接触刚度kτ以及界面摩擦系数μ等。同时划分成满足计算精度要求的网络。
S2、以网格中心为粒子坐标，网格体积为粒子体积V，生成刚性弹体粒子离散模型以及土体粒子离散模型。
赋值质量m＝ρV。设定土体粒子集合Bs和刚性弹体离散集合Bm，将步骤S1中刚性弹体和土体变量输入至对应粒子离散模型中。
S3、初始化刚性弹体及土体的信息，分别计算土体粒子的运动方程以及刚性弹体粒子的运动方程。
设定刚性弹体侵彻土体姿态、六个刚体自由度(三个方向速度、三个方向旋转速度)以及土体的地应力分布。并通过SPH离散，建立土体粒子的运动方程以及刚性弹体粒子的运动方程。
具体的，土体粒子的运动方程计算方法包括：
依据张量理论中的符号规定，土体的变形控制方程形式，每个土体粒子a的离散质量守恒和动量守恒方程如下：
式中，表示为时间微分；ma、ρa分别表示土体粒子a的质量和质量密度；va表示土体粒子a的速度矢量；/>表示侵彻外力；/>表示内部约束力；b表示土体粒子b；Ωa表示土体粒子a的支持域；Vb表示土体粒子b的体积；vb表示土体粒子b的速度矢量。gab土体粒子a对土体粒子b表示辅助核函数导数的离散形式，具体公式表示为gab＝g(xa-xb,ls)；g表示为辅助核函数梯度；xa表示粒子a的坐标；xb表示粒子b的坐标；ls表示为光滑长度。
其中，内部约束力取决于支持域内各点处变形梯度及本构关系，体现粒子与其支持域内其他粒子的相互作用，如下式：
式中，Va表示土体粒子a的体积；σa表示土体粒子a的柯西应力；σb表示土体粒子b的柯西应力；gba土体粒子b对土体粒子a表示辅助核函数导数的离散形式，具体公式表示为gba＝g(xb-xa,ls)；g表示为辅助核函数梯度；xa表示粒子a的坐标；xb表示粒子b的坐标；ls表示为光滑长度。
由于SPH法计算结果常表现出不稳定的数值震荡，因此可以采用人工粘度增加内部约束，仅对应力进行修正：
式中，αI表示人工粘度参数；cart表示人工声速；ls表示光滑长度；hac表示土体粒子a对土体粒子c的辅助核函数的离散形式，hbc表示土体粒子b对土体粒子c的辅助核函数的离散形式，I为二阶单位张量。c表示土体粒子c，vc表示土体粒子c的速度，Vc表示土体粒子c的体积。
为了能够得到更高精度的预测结果，SPH核函数hab＝h(xa-xb,l)及其导数gab＝g(xa-xb,l)的离散形式则根据周围粒子的坐标信息(xba＝[xb-xa,yb-ya,zb-za]T，式中，xa、xb分别表示土体粒子a，b的x轴坐标，ya、yb分别表示土体粒子a，b的y轴坐标，za、zb分别表示土体粒子a，b的z轴坐标。)求解得到，具体如下式：
式中，gab＝[gx,ab gy,ab gz,ab]T分别表示辅助核函数导数的x，y，z轴方向的离散形式；Wab，为原始核函数及三个方向导数的离散形式，并采用高斯函数作为SPH核函数，具体表达式为：
土体在被刚性弹体侵彻过程中发生形变，需要考虑土体强度对抵抗刚性弹体侵彻的能力，因此需要确定土体的本构关系，即，考虑土体材料参数与柯西应力σa大小。本实施例中，采用Drucker-Prager本构模型，引入土体材料参数内摩擦角粘聚力cm，剪胀角ψ对于应力值的影响。
假设土体的本构关系符合Drucker-Prager本构模型，根据客观应变力关系以及弹塑性理论，本构关系表达式如下：
式中，为柯西应力的客观应变率，D表示四阶弹性张量，m表示预定义的塑性势函数应力梯度，/>表示弹塑性理论中的塑性乘子，/>表示梯度算子，/>表示速度。
Drucker-Prager本构模型的屈服函数以及塑性势函数表达式为：
g＝q-tanψp，其中，p＝-tr(σ)/3；/>s＝σ-pI.
式中，tr(·)表示张量的迹。
塑性乘子从一致性条件中实例化，表达式为：
式中，K和G分别为体积模量和剪切模量，表示应变速率。
本实施例中，刚性弹体粒子的运动方程计算方法包括：
设定刚性弹体重心为x0＝(x0,y0,z0)，其六个自由度包含三个方向的速度v，三个方向的角速度ω0，刚性弹体质量为M0，三个方向的惯性矩为J0，刚性弹体重力为G，Text表示抗侵彻力；Sext表示刚性弹体受到的抗侵彻外合力矩；根据牛顿第二定律可得到合理与合力矩力学方程，表示式如下：
式中，v0表示刚性弹体的速度，ω0表示刚性弹体的角速度。
通过数值积分后，刚性弹体粒子d所对应的速度状态信息可通过刚性弹体的速度和角速度获得，公式为：
vd＝v0+ω0×(xd-x0)
式中，ω0表示刚性弹体的角速度。xd为刚性弹体粒子d的坐标，x0为刚性弹体重心坐标。
S4、计算刚性弹体与土体之间的接触摩擦力，基于接触摩擦力、土体粒子的运动方程以及刚性弹体粒子的运动方程计算刚性弹体的运动轨迹，记录刚性弹体速度为0时的侵彻水平深度和竖直深度。
如图2所示，为土体粒子与刚性弹体粒子之间的接触关系示意图。接触摩擦力包括：法向接触力fn和切向接触摩擦力fτ，并分别作用累加到土体粒子受到的侵彻外力、以及刚性弹体受到的抗侵彻力及抗侵彻外合力矩。
本实施例中，采用点对体积的界面接触摩擦模型，根据土体离散的土体粒子对于刚性弹体的覆盖体积计算得到接触摩擦力。
计算刚性弹体与土体接触界面受力大小，首先需要对于每个刚性弹体附近的土体粒子a∈Bs，利用其支持域Ωa内所包含的刚性弹体粒子b∈Ωa∩Bm，定义土体粒子嵌入刚性弹体的距离表达式：
式中，w0为调整土体粒子嵌入刚性弹体的嵌入量，αls为支持域半径。wa∈[0,1]是根据土体粒子支持域内的刚性弹体粒子积分，用以计算土体粒子与刚性弹体之间的距离，表达式为：
根据土体粒子对于刚性弹体的嵌入量则可以得到土体粒子对于刚性弹体的法向接触力，表达式为：
式中，kn为法向接触刚度；代表接触面外法向单位矢量，根据土体粒子支持域内的刚性弹体粒子，其表达式为：
式中，表示原核函数导数的三个方向。
土体粒子对于刚性弹体的切向接触力fτ,a为累积量，根据计算过程的累积，表达式为：
式中，kτ为切向接触刚度，Δt表示时间步，vτ,a为土体粒子相对刚性弹体的切向速度，公式为：
当发生滑动时，根据界面摩擦系数μ，计算得到切向接触摩擦力：
式中，为单位切向力方向。
根据上述接触力计算公式，可得土体粒子受到的侵彻外力：
式中，fn,a表示法向接触力；fτ,a表示切向接触摩擦力；
根据作用力和反作用力，得到刚性弹体受到的抗侵彻力：
式中，a表示土体粒子；Bs表示土体粒子集合；
刚性弹体受到的抗侵彻外合力矩为：
式中，xa0表示刚性弹体重心到土体粒子a的距离。
根据上述土体粒子的运动方程、刚性弹体粒子的运动方程以及设置的初始条件，计算，当刚性弹体侵彻土体过程静止，则认为侵彻完成。
实施例二：
本实施例提供一种刚性弹体侵彻多层土体深度的数值仿真预测系统，包括：划分单元、生成单元、计算单元及轨迹单元。
划分单元用于构建刚性弹体三维几何模型以及土体分布几何模型，并对刚性弹体三维几何模型以及土体分布几何模型进行参数赋值，以及网格划分。
根据实际情况创建不同分层形状的土体分布几何模型以及刚性弹体三维几何模型；对每一层土体分别赋值不同密度ρ，以及其材料本构参数：弹性模量E、泊松比ν、内摩擦角黏聚力cm以及剪胀角ψ等。设置刚性弹体与土体之间的法向接触刚度kn、切向接触刚度kτ以及界面摩擦系数μ等。同时划分成满足计算精度要求的网络。
生成单元用于以网格中心为粒子坐标，网格体积为粒子体积，生成刚性弹体粒子离散模型以及土体粒子离散模型。
赋值质量m＝ρV。设定土体粒子集合Bs和刚性弹体离散集合Bm，将步骤S1中刚性弹体和土体变量输入至对应粒子离散模型中。
计算单元用于初始化刚性弹体及土体的信息，分别计算土体粒子的运动方程以及刚性弹体粒子的运动方程。
设定刚性弹体侵彻土体姿态、六个刚体自由度(三个方向速度、三个方向旋转速度)以及土体的地应力分布。并通过SPH离散，建立土体粒子的运动方程以及刚性弹体粒子的运动方程。
具体的，土体粒子的运动方程计算方法包括：
依据张量理论中的符号规定，土体的变形控制方程形式，每个土体粒子a的离散质量守恒和动量守恒方程如下：
式中，表示为时间微分；ma、ρa分别表示土体粒子a的质量和质量密度；va表示土体粒子a的速度矢量；/>表示侵彻外力；/>表示内部约束力；b表示土体粒子b；Ωa表示土体粒子a的支持域；Vb表示土体粒子b的体积；vb表示土体粒子b的速度矢量。gab表示土体粒子a对土体粒子b表示辅助核函数导数的离散形式，具体公式表示为gab＝g(xa-xb,ls)；g表示为辅助核函数梯度；xa表示粒子a的坐标；xb表示粒子b的坐标；ls表示为光滑长度。
其中，内部约束力取决于支持域内各点处变形梯度及本构关系，体现粒子与其支持域内其他粒子的相互作用，如下式：
式中，Va表示土体粒子a的体积；σa表示土体粒子a的当前状态柯西应力；σb表示土体粒子b的当前状态柯西应力；gba表示土体粒子b对土体粒子a的辅助核函数导数的离散形式，具体公式表示为gba＝g(xb-xa,ls)；g表示为辅助核函数梯度；xa表示粒子a的坐标；xb表示粒子b的坐标；ls表示为光滑长度。
由于SPH法计算结果常表现出不稳定的数值震荡，因此可以采用人工粘度增加内部约束，仅对应力进行修正：
式中，αI表示人工粘度参数；cart表示人工声速；ls表示光滑长度；hac表示光滑粒子流体动力学法核函数的离散形式，hbc表示土体粒子b对土体粒子c的辅助核函数离散形式，I为二阶单位张量。c表示土体粒子c，vc表示土体粒子c的速度，Vc表示土体粒子c的体积。
为了能够得到更高精度的预测结果，SPH核函数hab＝h(xa-xb,l)及其导数gab＝g(xa-xb,l)的离散形式则根据周围粒子的坐标信息(xba＝[xb-xa,yb-ya,zb-za]T，式中，xa、xb分别表示土体粒子a，b的x轴坐标，ya、yb分别表示土体粒子a，b的y轴坐标，za、zb分别表示土体粒子a，b的z轴坐标。)求解得到，具体如下式：
式中，gab＝[gx,ab gy,ab gz,ab]T分别表示辅助核函数导数的x，y，z轴方向的离散形式；Wab，为原始核函数及三个方向导数的离散形式，并采用高斯函数作为SPH核函数，具体表达式为：
土体在被刚性弹体侵彻过程中发生形变，需要考虑土体强度对抵抗刚性弹体侵彻的能力，因此需要确定土体的本构关系，即，考虑土体材料参数与柯西应力σa大小。本实施例中，采用Drucker-Prager本构模型，引入土体材料参数内摩擦角粘聚力c，剪胀角ψ对于应力值的影响。
假设土体的本构关系符合Drucker-Prager本构模型，根据客观应变力关系以及弹塑性理论，本构关系表达式如下：
式中，为柯西应力的客观应变率，D表示四阶弹性张量，m表示预定义的塑性势函数应力梯度，/>表示弹塑性理论中的塑性乘子，/>表示梯度算子，/>表示速度。
Drucker-Prager本构模型的屈服函数以及塑性势函数表达式为：
g＝q-tanψp，
其中，p＝-tr(σ)/3；s＝σ-pI.
式中，tr()表示张量的迹。
塑性乘子从一致性条件中实例化，表达式为：
式中，K和G分别为体积模量和剪切模量，表示应变速率。
本实施例中，刚性弹体粒子的运动方程计算方法包括：
设定刚性弹体重心为x0＝(x0,y0,z0)，其六个自由度包含三个方向的速度v，三个方向的角速度ω0，刚性弹体质量为M0，三个方向的惯性矩为J0，刚性弹体重力为G，Text表示抗侵彻力；Sext表示刚性弹体受到的抗侵彻外合力矩；根据牛顿第二定律可得到合理与合力矩力学方程，表示式如下：
式中，v0表示刚性弹体的速度，ω0表示刚性弹体的角速度。
通过数值积分后，刚性弹体粒子d所对应的速度状态信息为：
vd＝v0+ω0×(xd-x0)，
式中，ω0表示刚性弹体的角速度。xd为刚性弹体粒子d的坐标，x0为刚性弹体重心坐标。
轨迹单元用于计算刚性弹体与土体之间的接触摩擦力，基于接触摩擦力、土体粒子的运动方程以及刚性弹体粒子的运动方程计算刚性弹体的运动轨迹，记录刚性弹体速度为0时的侵彻水平深度和竖直深度。
接触摩擦力包括：所述接触摩擦力包括：法向接触力fn和切向接触摩擦力fτ，并分别作用累加到土体粒子受到的侵彻外力、以及刚性弹体受到的抗侵彻力及抗侵彻外合力矩。
本实施例中，采用点对体积的界面接触摩擦模型，根据土体离散的土体粒子对于刚性弹体的覆盖体积计算得到接触摩擦力。
计算刚性弹体与土体接触界面受力大小，首先需要对于每个刚性弹体附近的土体粒子a∈Bs，利用其支持域Ωa内所包含的刚性弹体粒子b∈Ωa∩Bm，定义土体粒子嵌入刚性弹体的距离表达式：
式中，w0为调整土体粒子嵌入刚性弹体的嵌入量，αls为支持域半径。wa∈[0,1]是根据土体粒子支持域内的刚性弹体粒子积分，用以计算土体粒子与刚性弹体之间的距离，表达式为：
根据土体粒子对于刚性弹体的嵌入量则可以得到土体粒子对于刚性弹体的法向接触力，表达式为：
式中，kn为法向接触刚度；代表接触面外法向单位矢量，根据土体粒子支持域内的刚性弹体粒子，其表达式为：
式中，表示原核函数导数的三个方向。
土体粒子对于刚性弹体的切向接触力fτ,a为累积量，根据计算过程的累积，表达式为：
fτ,a＝∫Δtkτvτ,adt，
式中，kτ为切向接触刚度，Δt表示时间增量步，vτ,a为土体粒子相对刚性弹体的切向速度，公式为：
当发生滑动时，根据界面摩擦系数μ，计算得到切向接触摩擦力：
式中，为单位切向力方向。
根据上述接触力计算公式，可得土体粒子受到的侵彻外力：
式中，fn,a表示法向接触力；fτ,a表示切向接触摩擦力；
根据作用力和反作用力，得到刚性弹体受到的抗侵彻力：
式中，a表示土体粒子；Bs表示土体粒子集合；
刚性弹体受到的抗侵彻外合力矩为：
式中，xa0表示刚性弹体重心到土体粒子a的距离。
根据上述土体粒子的运动方程、刚性弹体粒子的运动方程以及设置的初始条件，计算，当刚性弹体侵彻土体过程静止，则认为侵彻完成。
实施例三：
为了验证本发明提出的方法的效果，本实施例提供一种测试算例。
如图3所示，刚性弹体为圆柱体，直径180mm，长度800mm，重量共75kg，侵彻姿态γ＝50°，侵彻竖向速度为100km/h，水平速度250km/h，如表1所示，为调整前两层土体材料参数。图4为初始参数下刚性弹体侵彻的最终侵彻深度结果。
表1
若有试验结果对比，可根据调整材料参数、接触摩擦系数进行再一次预测，直至满足误差要求为止。经过刚性弹体和土体参数条件，调整参数如下，刚性弹体直径190mm，长度1200mm，重量共85kg，侵彻姿态γ＝55°，侵彻竖向速度为200km/h，水平速度300km/h，如表所示材料参数。调整后参数如表2所示，其预测结果如图5所示。
表2
本发明所提出的方法模拟了刚性弹体侵彻土体的全部过程，综合考虑土体强度、刚性弹体土体界面接触以及三维刚性刚性弹体姿态和6个自由度等主要影响侵彻深度的因素，能够实现刚性弹体侵彻土体的深度分析，能够有效预测侵彻的深度，进而能够指导刚性弹体和土体接触界面侵彻的抑制以及参数的选择。
以上所述的实施例仅是对本发明优选方式进行的描述，并非对本发明的范围进行限定，在不脱离本发明设计精神的前提下，本领域普通技术人员对本发明的技术方案做出的各种变形和改进，均应落入本发明权利要求书确定的保护范围内。

Claims (4)

1.一种刚性弹体侵彻多层土体深度的数值仿真预测方法，其特征在于，包括以下步骤：
S1、构建刚性弹体三维几何模型以及土体分布几何模型，并对所述刚性弹体三维几何模型以及土体分布几何模型进行参数赋值，以及网格划分；
S2、以网格中心为粒子坐标，网格体积为粒子体积，生成刚性弹体粒子离散模型以及土体粒子离散模型；
S3、初始化刚性弹体及土体的信息，分别计算土体粒子的运动方程以及刚性弹体粒子的运动方程；
S4、计算刚性弹体与土体之间的接触摩擦力，基于所述接触摩擦力、所述土体粒子的运动方程以及所述刚性弹体粒子的运动方程计算刚性弹体的运动轨迹，记录刚性弹体速度为0时的侵彻水平深度和竖直深度；
所述土体粒子的运动方程包括：
，
式中，表示为时间微分；m a 、ρ a 分别表示土体粒子a的质量和质量密度；v a 表示土体粒子a的速度矢量；/>表示侵彻外力；/>表示内部约束力；/>土体粒子a对土体粒子b表示辅助核函数导数的离散形式；b表示土体粒子b；Ω a 表示土体粒子a的支持域；V b 表示土体粒子b的体积；v b 表示土体粒子b的速度矢量；
所述刚性弹体粒子的运动方程包括：
，
式中，表示为时间微分；M 0表示刚性弹体粒子质量；v 0表示刚性弹体的速度；T ext表示刚性弹体受到的抗侵彻力；G表示刚性弹体粒子重力；J 0表示刚性弹体三个方向上的惯性矩；S ext表示刚性弹体受到的抗侵彻外合力矩；ω 0表示刚性弹体的三个方向上的角速度。
2.根据权利要求1所述刚性弹体侵彻多层土体深度的数值仿真预测方法，其特征在于，所述接触摩擦力包括：法向接触力f n 和切向接触摩擦力f τ ，并分别作用累加到土体粒子受到的侵彻外力、以及刚性弹体受到的抗侵彻力及抗侵彻外合力矩；
所述土体粒子受到的侵彻外力包括：
，
式中，f n,a 表示法向接触力；f τ,a 表示切向接触摩擦力；
所述刚性弹体受到的抗侵彻力包括：
，
式中，a表示土体粒子；B s 表示土体粒子集合；
所述刚性弹体受到的抗侵彻外合力矩包括：
，
式中，x a0表示刚性弹体重心到土体粒子a的距离。
3.一种刚性弹体侵彻多层土体深度的数值仿真预测系统，所述系统用于实现如权利要求1-2任一项所述的方法，其特征在于，包括：划分单元、生成单元、计算单元及轨迹单元；
所述划分单元用于构建刚性弹体三维几何模型以及土体分布几何模型，并对所述刚性弹体三维几何模型以及土体分布几何模型进行参数赋值，以及网格划分；
所述生成单元用于以网格中心为粒子坐标，网格体积为粒子体积，生成刚性弹体粒子离散模型以及土体粒子离散模型；
所述计算单元用于初始化刚性弹体及土体的信息，分别计算土体粒子的运动方程以及刚性弹体粒子的运动方程；
所述轨迹单元用于计算刚性弹体与土体之间的接触摩擦力，基于所述接触摩擦力、所述土体粒子的运动方程以及所述刚性弹体粒子的运动方程计算刚性弹体的运动轨迹，记录刚性弹体速度为0时的侵彻水平深度和竖直深度；
所述土体粒子的运动方程包括：
，
式中，表示为时间微分；m a 、ρ a 分别表示土体粒子a的质量和质量密度；v a 表示土体粒子a的速度矢量；/>表示侵彻外力；/>表示内部约束力；/>土体粒子a对土体粒子b表示辅助核函数导数的离散形式；b表示土体粒子b；Ω a 表示土体粒子a的支持域；V b 表示土体粒子b的体积；v b 表示土体粒子b的速度矢量；
所述刚性弹体粒子的运动方程包括：
，
式中，表示为时间微分；M 0表示刚性弹体质量；v 0表示刚性弹体的速度；T ext表示刚性弹体受到的抗侵彻力；G表示刚性弹体重力；J 0表示刚性弹体三个方向上的惯性矩；S ext表示刚性弹体受到的抗侵彻外合力矩；ω 0表示刚性弹体的三个方向上的角速度。
4.根据权利要求3所述刚性弹体侵彻多层土体深度的数值仿真预测系统，其特征在于，所述接触摩擦力包括：法向接触力f n 和切向接触摩擦力f τ ，并分别作用累加到土体粒子受到的侵彻外力、以及刚性弹体受到的抗侵彻力及抗侵彻外合力矩；
所述土体粒子受到的侵彻外力包括：
，
式中，f n,a 表示法向接触力；f τ,a 表示切向接触摩擦力；
所述刚性弹体受到的抗侵彻力包括：
，
式中，a表示土体粒子；B s 表示土体粒子集合；
所述刚性弹体受到的抗侵彻外合力矩包括：
，
式中，x a0表示刚性弹体重心到土体粒子a的距离。
CN202311609709.7A 2023-11-28 2023-11-28 刚性弹体侵彻多层土体深度的数值仿真预测方法及系统 Active CN117610393B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202311609709.7A CN117610393B (zh) 2023-11-28 2023-11-28 刚性弹体侵彻多层土体深度的数值仿真预测方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202311609709.7A CN117610393B (zh) 2023-11-28 2023-11-28 刚性弹体侵彻多层土体深度的数值仿真预测方法及系统

Publications (2)

Publication Number Publication Date
CN117610393A CN117610393A (zh) 2024-02-27
CN117610393B true CN117610393B (zh) 2024-05-14

Family

ID=89943869

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202311609709.7A Active CN117610393B (zh) 2023-11-28 2023-11-28 刚性弹体侵彻多层土体深度的数值仿真预测方法及系统

Country Status (1)

Country Link
CN (1) CN117610393B (zh)

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111666678A (zh) * 2020-06-03 2020-09-15 西安近代化学研究所 一种半流体侵彻的二维建模方法
CN112084659A (zh) * 2020-09-09 2020-12-15 北京理工大学 考虑侵蚀效应的弹体高速冲击混凝土侵彻性能的预测方法
CN115422750A (zh) * 2022-09-02 2022-12-02 山东非金属材料研究所 一种获取刚性弹体弹速与侵彻深度关系的试验方法
CN116580792A (zh) * 2023-05-10 2023-08-11 中国人民解放军火箭军工程设计研究院 聚能装药侵彻深度预估方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111666678A (zh) * 2020-06-03 2020-09-15 西安近代化学研究所 一种半流体侵彻的二维建模方法
CN112084659A (zh) * 2020-09-09 2020-12-15 北京理工大学 考虑侵蚀效应的弹体高速冲击混凝土侵彻性能的预测方法
CN115422750A (zh) * 2022-09-02 2022-12-02 山东非金属材料研究所 一种获取刚性弹体弹速与侵彻深度关系的试验方法
CN116580792A (zh) * 2023-05-10 2023-08-11 中国人民解放军火箭军工程设计研究院 聚能装药侵彻深度预估方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
光滑粒子流体动力学(SPH)中摩擦接触的模拟;王建;吴浩;顾冲时;花卉;;中国科学:技术科学;20131120(第11期);第1209-1218页 *
刚性弹体对素混凝土厚靶侵彻响应的LDPM数值模拟研究;冯君;李文彬;徐磊;陈宇;陈曦;;振动与冲击;20180428(第08期);正文全文 *
基于FEM/SPH算法的弹丸侵彻土体数值模拟;罗伟铭;石少卿;田镇华;尹平;;后勤工程学院学报;20140330(第02期);第13-17页 *
弹丸侵彻混凝土的SPH算法;宋顺成, 才鸿年;爆炸与冲击;20030225(第01期);第56-60页 *

Also Published As

Publication number Publication date
CN117610393A (zh) 2024-02-27

Similar Documents

Publication Publication Date Title
Henke Influence of pile installation on adjacent structures
Yu et al. Multiscale method for long tunnels subjected to seismic loading
CN109241588A (zh) 一种基于拟连续地质力学模型的单裂缝扩展的模拟方法
Gharehdash et al. Numerical modeling of the dynamic behaviour of tunnel lining in shield tunneling
CN108572401B (zh) 缝洞组合模型的构建方法及探测储层缝洞变形的方法
Marina et al. Simulation of the hydraulic fracturing process of fractured rocks by the discrete element method
Chen et al. Numerical study of 3-D liquid sloshing in an elastic tank by MPS-FEM coupled method
CN104239637A (zh) 一种离散元爆堆形态模拟方法
Li et al. A new dynamic modelling methodology of a force measuring system for hypersonic impulse wind tunnel
Liang et al. The imposition of nonconforming Neumann boundary condition in the material point method without boundary representation
CN113033060B (zh) 一种用于预测复杂煤层开采结构的优化方法
Kawamoto et al. A review of numerical analysis of tunnels in discontinuous rock masses
CN117610393B (zh) 刚性弹体侵彻多层土体深度的数值仿真预测方法及系统
Wobbes et al. Modeling of liquefaction using two-phase FEM with UBC3D-PLM model
Xiong et al. Development of an unresolved CFD-DEM method for interaction simulations between large particles and fluids
Valdi et al. Numerical investigation of water entry problem of pounders with different geometric shapes and drop heights for dynamic compaction of the seabed
CN103206203A (zh) 油井单一射孔出砂的分析方法
Santos et al. Compression and shear-wave velocities in discrete particle simulations of quartz granular packings: Improved Hertz-Mindlin contact model
Heymsfield et al. Assessment of soil modeling capability for orion contingency land landing
CN114722681A (zh) 一种盾构施工引发地面沉降模拟预测方法
Don et al. A digital twinning methodology for vibration prediction and fatigue life prognosis of vertical oil well Drillstrings
CN102493800B (zh) 一种射孔弹性能参数的欧拉获取方法
Zhang et al. Rotation errors in numerical manifold method and a correction based on large deformation theory
Te Kamp et al. Micromechanical back-analysis of laboratory tests on rock
Mohammad et al. Erosion and Structural Integrity of Mud Pulse Telemetry Tools: Numerical Simulation and Field Studies

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant