CN113032974A - 一种六面体单元有限质点法 - Google Patents

一种六面体单元有限质点法 Download PDF

Info

Publication number
CN113032974A
CN113032974A CN202110246003.3A CN202110246003A CN113032974A CN 113032974 A CN113032974 A CN 113032974A CN 202110246003 A CN202110246003 A CN 202110246003A CN 113032974 A CN113032974 A CN 113032974A
Authority
CN
China
Prior art keywords
particle
unit
mass
hexahedral
displacement
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.)
Pending
Application number
CN202110246003.3A
Other languages
English (en)
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.)
Tsinghua University
Original Assignee
Tsinghua University
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 Tsinghua University filed Critical Tsinghua University
Priority to CN202110246003.3A priority Critical patent/CN113032974A/zh
Publication of CN113032974A publication Critical patent/CN113032974A/zh
Pending legal-status Critical Current

Links

Images

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
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/17—Mechanical parametric or variational design
    • 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)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • Databases & Information Systems (AREA)
  • Computer Hardware Design (AREA)
  • Algebra (AREA)
  • Evolutionary Computation (AREA)
  • Software Systems (AREA)
  • Computing Systems (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Operations Research (AREA)
  • Probability & Statistics with Applications (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

一种六面体单元有限质点法，包括以下步骤：首先对模型参数初始化，并更新质点位移，然后对六面体单元逆向运动求解质点纯变形，进而计算出质点内力，最终带入中央差分运动公式求解下一步质点位移。本发明通过对传统有限质点法的六面体单元逆向转动以及建立局部坐标系步骤进行改进，采用迭代计算对质点的位移和内力进行更新，使其进一步提高了计算准确度，并在保证计算精度的前提下提高计算效率，有效降低了刚体转动高阶误差，简化了计算流程，也使得推导过程更加简单清晰。本发明的推导过程与方法同样适用于其他类型的实体单元，具有普遍适用意义。

Description

一种六面体单元有限质点法
技术领域
本发明属于结构力学领域，具体的涉及一种六面体单元有限质点法。
背景技术
有限质点法是基于质点运动和向量式结构力学提出的面向结构分析的数值分析方法，对于结构的不连续性以及非线性问题，也能进行准确分析；该方法不再使用连续体的概念分析对象，而是基于牛顿运动定律，将真实的物理模型离散成具有有限质点数的质点系，分析质点之间相互独立的运动；将结构运动变形时间分解为有限个时间段，每一个时间段就是一个途径单元，途径单元长度为迭代步长，在每一个途径单元内对各个质点的受力和运动情况计算，迭代计算质点位置，求得整体结构响应。有限质点法采用逆向运动消除结构刚体位移对质点纯变形计算的影响，对线性、强非线性和不连续变形问题的分析都是基于统一的框架，无需改变计算流程或者控制方程，并且在计算过程中无需进行刚度矩阵的集成以及求解复杂的非线性方程组。
对于结构分析，通过建立实体模型反映的响应更为准确，所需的实体单元类型主要有六面体和四面体单元，相对于四面体单元，六面体单元有着更高的计算精度，更广的应用范围。但是当前有限质点法的研究还处于初步阶段，对于杆、梁单元、平面膜单元、薄壳单元研究较多，对于实体单元，尤其是六面体单元的研究较少。当前传统六面体单元有限质点法中的逆向运动步骤是以六面体单元质点作为局部坐标系原点求解质点纯变形，虽然也消除大部分的单元刚体位移，残余的误差相对于质点纯变形也只是高阶误差，但是人为的将单元部分质点的纯变形归0，导致高阶误差分布不均，增加了计算误差，尽管传统六面体单元有限质点法通过单元的力矩平衡关系求解出了纯变形为0节点的质点内力，但是却增大了单元迭代计算过程中质点内力值的波动，需要更小的迭代步长来防止中央差分运动公式的发散，增加了迭代次数，降低了计算效率。并且当前传统六面体单元有限质点法在平面内转动步骤中多次判断转角正负和向量求模，使得计算过程较为复杂。
发明内容
本发明的目的在于提供一种新的六面体单元有限质点法，使其进一步提高计算准确度，并在保证计算精度的前提下提高计算效率，降低刚体转动高阶误差，有效简化计算流程。
本发明采用如下的技术方案：
一种六面体单元有限质点法，包括如下步骤：
步骤1、将六面体单元质量平均分配到8个质点上，每个质点的质量为m，质点与质点之间不具备质量，以任意一点为全局坐标系O-xyz原点，坐标轴单位向量为：ex＝[1 0 0]T，ey＝[0 1 0]T，ez＝[0 0 1]T，质点的初速度为
Figure BDA0002964108430000025
初位移为x0、初始内力为f0和初始外力为F0，对分析时长ET进行分段，每一个时间段是一个分析步长h，所需要的总步数为N；
步骤2、计算第-1步的虚拟质点位移x-1：
Figure BDA0002964108430000021
式中，x-1为第-1步的虚拟质点位移；ζ为阻尼系数；P0＝f0+F0为初始质点外力和；
步骤3、根据上步计算结果，更新质点位移：将第n步和第n-1步的质点位移xn和xn-1作为计算质点位移xn+1的参数；
步骤4、对第n+1步的六面体单元逆向运动得到质点的纯变形，包括以下步骤：
步骤4-1、对第n+1步的六面体单元进行逆向平移：以第n+1步的六面体单元质心为参考点，整体逆向平移至第n步的六面体单元质心上；
步骤4-2、对步骤4-1中经过逆向平移的六面体单元进行第一次逆向转动：将经过逆向平移的六面体单元水平方向上的四个平面质心构造参考平面求解法向量，同理，对第n步的六面体单元也构造出参考平面法向量，以两个法向量的垂直向量为转动向量，夹角为转动角，经过逆向平移的六面体单元质心为转动原点进行整体逆向转动；
步骤4-3、对步骤4-2中经过第一次逆向转动的六面体单元进行第二次逆向转动：将经过第一次逆向转动的六面体单元垂直方向上的四个平面质心构造参考平面求解法向量，同理，对第n步的六面体单元也构造出参考平面法向量，以两个法向量的垂直向量为转动向量，夹角为逆向转动角，经过第一次逆向转动的六面体单元质心为转动原点进行整体逆向转动；
将经过第二次逆向转动的六面体单元质点位置与第n步的六面体单元质点位置作差，计算出质点纯变形：
Figure BDA0002964108430000022
式中，qi为六面体单元各质点的纯变形；
Figure BDA0002964108430000023
为第n步的各质点位置；xi″′为经过逆向运动后的各质点位置；
步骤5、计算质点内力：以第n步六面体单元质心为原点建立局部坐标系
Figure BDA0002964108430000024
将全局坐标系O-xyz下各个质点位置和纯变形转换到局部坐标系中，利用单元形函数求解出虚拟质点内力，将虚拟质点内力正向运动，获得实际位置的质点内力；
步骤6、根据质点中央差分运动公式计算质点位移xn+1：
Figure BDA0002964108430000031
式中，Pn＝fn+Fn为第n步质点外力和，fn为第n步质点内力，Fn为第n步质点外力；
步骤7、判断分析步数是否达到总步数，当n+1＝N时，达到所需总步数，循环计算结束，否则转达步骤3。
本发明与传统有限质点法相比，其显著优点为：1)本发明使用六面体单元各平面质心建立逆向转动所需的参考平面，加强参考平面与六面体单元各质点的联系，提高计算准确度；2)本发明采用的逆向运动方式有效降低了节点内力值的波动性，在保证计算精度的前提下选取更大的迭代步长，提高了计算效率；3)本发明使用单元质心作为局部坐标系原点，使六面体单元各质点纯变形分布更均匀，降低了刚体转动高阶误差；4)本发明使用质点中央差分运动公式求解质点位移，有效避免了有限元法计算质点内力过程中的集成整体刚度矩阵步骤，简化了计算流程。
附图说明
图1是一种六面体单元有限质点法流程图。
图2是六面体单元运动和变形示意图。
图3是逆向平移后的六面体单元位置示意图。
图4是第一次逆向转动后的六面体单元位置示意图。
图5是第二次逆向转动后的六面体单元位置示意图。
图6是八节点六面体单元形函数坐标系示意图。
图7是悬臂梁模型受力分析结果图。
具体实施方式
下面结合附图对本发明作进一步详细描述。
本发明提供的一种六面体单元有限质点法，其包括以下步骤：
步骤1、将六面体单元质量平均分配到8个质点上，每个质点的质量为m，质点与质点之间不具备质量，以任意一点为全局坐标系O-xyz原点，坐标轴单位向量为：ex＝[1 0 0]T，ey＝[0 1 0]T，ez＝[0 0 1]T，质点的初速度为
Figure BDA0002964108430000041
初位移为x0、初始内力为f0和初始外力为F0，对分析时长ET进行分段，每一个时间段是一个分析步长h，所需要的总步数为N；
步骤2、计算第-1步的虚拟质点位移x-1：
Figure BDA0002964108430000042
式中，x-1为第-1步的虚拟质点位移；ζ为阻尼系数；P0＝f0+F0为初始质点外力和；
步骤3、根据上步计算结果，更新质点位移：将第n步和第n-1步的质点位移xn和xn-1作为计算第n+1步质点位移xn+1的参数；
步骤4、对第n+1步的六面体单元逆向运动得到质点的纯变形：首先对六面体单元进行逆向平移，再进行第一次逆向转动，最后进行第二次逆向转动，并与第n步质点位移比较，计算出质点纯变形：
Figure BDA0002964108430000043
式中，qi为六面体单元各质点的纯变形；
Figure BDA0002964108430000044
为第n步的各质点位置；xi″′为经过逆向运动后的各质点位置；
步骤5、计算质点内力：以第n步六面体单元质心为原点建立局部坐标系
Figure BDA0002964108430000045
将全局坐标系O-xyz下各个质点位置和纯变形转换到局部坐标系中，利用单元形函数求解出虚拟质点内力，将虚拟质点内力正向运动，获得实际位置的质点内力；
步骤6、根据质点中央差分运动公式计算质点位移xn+1：
Figure BDA0002964108430000046
式中，Pn＝fn+Fn为第n步质点外力和，fn为第n步质点内力，Fn为第n步质点外力；
步骤7、判断分析步数是否达到总步数，当n+1＝N时，达到所需总步数，循环计算结束，否则转达步骤3。
本发明有效降低了计算质点纯变形的误差，并且在保证计算精度的前提下提高计算效率。
下面进行更详细的描述。
参见图1，一种六面体单元有限质点法具体计算步骤如下：
第一步，将六面体单元离散为8个质点，并将六面体单元的质量平均分配至单元8个节点处的质点上，质点与质点之间不具备质量，质点质量的计算公式为：
Figure BDA0002964108430000051
式中，mi为质点i质量；M为六面体单元的质量；i＝a,b,c...g,h；
以任意一点为全局坐标系O-xyz原点，坐标轴单位向量为：ex＝[1 0 0]T，ey＝[0 10]T，ez＝[0 0 1]T，设定质点初位移
Figure BDA0002964108430000052
初速度
Figure BDA0002964108430000053
初始内力
Figure BDA0002964108430000054
和初始外力
Figure BDA0002964108430000055
i＝a,b,c...g,h，对分析时长ET进行分段，每一个时间段是一个途径单元，途径单元为tn到tn+1时刻，途径单元的长度是一个分析步长h，所需要的总步数为N；
接下来在一个途径单元tn到tn+1时刻内对六面体单元有限质点法步骤详细描述，对应的步数为第n到n+1步；
如图2所示，对第n到n+1步经过运动和变形的六面体单元进行如下定义：在全局坐标系下，第n步单元节点记为：an-bn-cn-dn-en-fn-gn-hn，第n+1步单元节点记为：an+1-bn+1-cn+1-dn+1-en+1-fn+1-gn+1-hn+1；
第二步，计算第-1步的虚拟质点位移
Figure BDA0002964108430000056
计算公式如下：
Figure BDA0002964108430000057
式中，
Figure BDA0002964108430000058
为第-1步的虚拟质点位移；ζ为阻尼系数；
Figure BDA0002964108430000059
为初始质点外力和；i＝a,b,c...g,h；
第三步，根据上步计算结果，更新质点位移：将第n步和第n-1步的质点位移
Figure BDA00029641084300000510
和
Figure BDA00029641084300000511
作为计算质点位移
Figure BDA00029641084300000512
的参数；
第四步，对第n+1步的六面体单元逆向运动得到质点的纯变形，逆向运动包括逆向平移、第一次逆向转动和第二次逆向转动。
逆向平移，如图3所示，将第n+1步的六面体单元以单元质心
Figure BDA00029641084300000513
为参考点，将六面体单元整体逆向平移，使得第n+1步单元质心
Figure BDA00029641084300000514
与第n步单元质心
Figure BDA00029641084300000515
重合，将经过逆向平移第n+1步的六面体单元节点记为：a′-b′-c′-d′-e′-f′-g′-h′，位置矢量x′i为：
Figure BDA0002964108430000061
第一次逆向转动，如图4所示，将经过逆向平移的六面体单元水平方向上四个平面的质心连接，构成两个三角形参考平面1’-4’-6’和2’-4’-6’，分别对两参考平面求法向量vec′1和vec′2，同理，根据第n步的六面体单元的参考平面1-4-6和2-4-6求解出参考平面法向量vec1和vec2；vec1与vec′1的夹角及其转动轴线向量分别为：
Figure BDA0002964108430000062
式中，θ1为vec1与vec′1的夹角；
Figure BDA0002964108430000063
式中，nvec1为vec1与vec′1的转动轴线向量；
同理，vec2与vec′2的夹角及其转动轴线向量分别为：θ2和nvec2；
对θ1与θ2、nvec1与nvec2分别取平均：
Figure BDA0002964108430000064
式中，θA为第一次逆向转动所需的转角；
Figure BDA0002964108430000065
式中，nvecA第一次逆向转动所需的转动轴线向量；
计算出第一次逆向转动的转动矩阵：
R′(-θΑ)＝[1-cos(-θΑ)]V′2+sin(-θΑ)V′ (式8)
式中，R′(-θΑ)为第一次逆向转动的转动矩阵；
Figure BDA0002964108430000066
Figure BDA0002964108430000067
以第n步的单元质心
Figure BDA0002964108430000068
为旋转中心对节点a′-b′-c′-d′-e′-f′-g′-h′逆向转动，获得新的单元节点位置为：a″-b″-c″-d″-e″-f″-g″-h″，位置矢量x″i为：
x″i＝R′(-θΑ)x′i(i＝a,b,c...g,h) (式9)
第二次逆向转动，如图5所示，将经过第一次逆向转动的六面体单元竖直方向上四个平面的质心连接，构成两个三角形参考平面5’-4’-6’和3’-4’-6’，分别对两参考平面求法向量vec′3和vec′4，同理，根据第n+1步的六面体单元的参考平面5-4-6和3-4-6求解出参考平面法向量vec3和vec4；vec3与vec′3的夹角及其转动轴线向量分别为：
Figure BDA0002964108430000071
式中，θ3为vec3与vec′3的夹角；
Figure BDA0002964108430000072
式中，nvec3为vec3与vec′3的转动轴线向量；
同理，vec4与vec′4的夹角及其转动轴线向量分别为：θ4和nvec4；
对θ3与θ4、nvec3与nvec4分别取平均：
Figure BDA0002964108430000073
式中，θB为第二次逆向转动所需的转角；
Figure BDA0002964108430000074
式中，nvecB第二次逆向转动所需的转动轴线向量；
计算出第二次逆向转动的转动矩阵：
R″(-θB)＝[1-cos(-θB)]V″2+sin(-θB)V″ (式14)
式中，R″(-θB)为第一次逆向转动的转动矩阵；
Figure BDA0002964108430000075
Figure BDA0002964108430000076
以tn时的单元质心
Figure BDA0002964108430000077
为旋转中心对节点a″-b″-c″-d″-e″-f″-g″-h″逆向转动，获得新的单元节点位置为：a″′-b″′-c″′-d″′-e″′-f″′-g″′-h″′，位置矢量x″′i为：
x″′i＝R″(-θB)xi″(i＝a,b,c...g,h) (式15)
将质点在第n步的位移矢量与经过逆向转动的质点位移矢量作差，得到六面体单元各质点的纯变形：
Figure BDA0002964108430000081
式中，qi为六面体单元各质点的纯变形；
需要特别说明的是，在整个逆向运动过程中，逆向平移是以六面体单元质心为参考点整体逆向平移，逆向转动中六面体单元的每一个质点都是绕同一旋转向量旋转相同的角度，保证了六面体单元的各个质点之间不会发生相对位移。
第五步、计算质点内力：如图5所示，以第n步单元质心作为原点建立局部坐标系
Figure BDA0002964108430000082
坐标轴方向与全局坐标系一致，将全局坐标系下质点位移矢量转换到局部坐标系中：
Figure BDA0002964108430000083
式中，
Figure BDA0002964108430000084
为在局部坐标系下的第n+1步质点位移矢量；
六面体单元质点内力计算将在局部坐标系下进行；
对第n步的六面体单元建立形函数，建立方式和表达形式和传统有限元方法相同，如图6所示，形函数坐标系为：O″-rst，八节点六面体单元的形函数表达式为：
Figure BDA0002964108430000085
式中，r,s,t分别为形函数表达式参数，ri＝±1，si＝±1，ti＝±1为形函数坐标系内质点i的坐标，i＝a,b,c...g,h；
将局部坐标系内坐标转换到形函数坐标系中，坐标系的转换通过雅各比矩阵J来实现：
Figure BDA0002964108430000086
式中，
Figure BDA0002964108430000087
为在局部坐标系内质点i的坐标，i＝a,b,c...g,h；
对雅克比矩阵求逆，得到形函数在局部坐标系中的导数为：
Figure BDA0002964108430000091
进而计算出六面体单元的应变增量
Figure BDA0002964108430000092
为：
Figure BDA0002964108430000093
式中，几何矩阵[B]＝[[B1][B2][B3][B4][B5][B6][B7][B8]]，
Figure BDA0002964108430000094
[q]＝[qa qb qc qd qe qf qgqh]；
六面体实体单元材料为弹性，本构关系矩阵D可由胡克定律得到：
Figure BDA0002964108430000095
式中，E为材料的弹性模量；μ为材料的泊松比；
Figure BDA0002964108430000096
六面体单元变形后的应力
Figure BDA0002964108430000097
为：
Figure BDA0002964108430000098
式中，
Figure BDA0002964108430000099
为第n步六面体单元的初始应力值；
六面体单元经过两次逆向转动的虚拟质点内力
Figure BDA00029641084300000910
为：
Figure BDA00029641084300000911
式中，虚拟质点内力
Figure BDA00029641084300000912
Figure BDA00029641084300000913
为虚拟质点i内力分量，i＝a,b,c...g,h；
将虚拟质点内力分量
Figure BDA00029641084300000914
正向运动，求得实际位置的质点内力分量
Figure BDA00029641084300000915
Figure BDA0002964108430000101
式中，R′(θA)和R″(θΒ)分别为第一次和第二次正向运动转动矩阵；
第六步、根据质点中央差分运动公式计算第n+1步质点i位移
Figure BDA0002964108430000102
Figure BDA0002964108430000103
式中，Pi n＝fi n+Fi n为六面体单元第n步质点i外力和，Fi n为六面体单元第n步质点i外力，i＝a,b,c...g,h；
第七步、判断分析步数是否达到总步数，当n+1＝N时，达到所需总步数，循环计算结束，否则转达步骤3。
下面结合实施例对本发明做进一步详细的描述。
实施例：
建立悬臂梁受力变形模型，模型参数及受力情况如表1所示：
表1悬臂梁参数及受力情况
Figure BDA0002964108430000104
对上述悬臂梁模型分别使用本发明所述的六面体单元有限质点法，有限元软件Abaqus，材料力学理论公式进行静力学分析，其中材料力学理论公式为：
Figure BDA0002964108430000111
式中，H为悬臂梁模型在竖直方向上产生的位移；F为悬臂梁模型自由末端受力；L为悬臂梁模型长度；E为悬臂梁模型材料弹性模量；
Figure BDA0002964108430000112
为转动惯量，b为悬臂梁模型截面的宽度，d为悬臂梁模型截面的高度；
表2各方法的计算结果
Figure BDA0002964108430000113
根据表2的计算结果看出本发明提出的六面体单元有限质点法与Abaqus软件和理论值的分析结果接近，验证了本发明提出的一种六面体单元有限质点法的可行性。
本实施例中采用的六面体单元有限质点法对悬臂梁模型建模使用了240个六面体单元，但是计算结果与使用960个六面体单元的Abaqus商业软件的计算结果十分接近，说明了本发明中的有限质点法对网格数量要求不高，只需对模型划分较少数量的网格即可得到精确的解。
本发明提出的六面体单元有限质点法，改进了传统有限质点法逆向运动和局部坐标系的设定步骤，其余部分还是遵循传统有限质点法的理论。对上述悬臂梁模型采用传统六面体单元有限质点法分析，因其对六面体单元质点纯变形的求解分布不均，导致所需迭代步长至少为1×10-6S，否则就会使质点内力值过大导致中央差分运动公式计算结果发散，并且传统六面体单元有限质点法在逆向运动和局部坐标系转换步骤计算更为复杂，因此采用传统六面体单元有限质点法对上述悬臂梁模型分析耗时至少是本发明提出的六面体单元有限质点法的10倍；
使用本发明中的有限质点法、Abaqus软件、材料力学理论公式，对悬臂梁模型受力后末端在竖直方向上位移变化结果如图7所示，有限质点法对静力学的分析，实际上是通过动态迭代计算分析直至收敛到静态解，而使用理论公式和Abaqus软件求解直接求出静力学分析结果，因此通过有限质点法对悬臂梁模型静力学分析结果是逐渐收敛的，这里的收敛速度没有实际的意义，通过调整阻尼系数来调节收敛速度。

Claims (3)

1.一种六面体单元有限质点法，其特征在于，包括以下步骤：
步骤1、将六面体单元质量平均分配到8个质点上，每个质点的质量为m，质点与质点之间不具备质量，以任意一点为原点建立全局坐标系O-xyz，坐标轴单位向量为：ex＝[1 0 0]T，ey＝[0 1 0]T，ez＝[0 0 1]T，质点的初速度为
Figure FDA0002964108420000011
初位移为x0、初始内力为f0和初始外力为F0，对分析时长ET进行分段，每一个时间段是一个分析步长h，所需要的总步数为N；
步骤2、计算第-1步的虚拟质点位移x-1：
Figure FDA0002964108420000012
式中，x-1为第-1步的虚拟质点位移；ζ为阻尼系数；P0＝f0+F0为初始质点外力和；
步骤3、根据上步的计算结果，更新质点位移；
步骤4、对第n+1步的六面体单元逆向运动得到质点的纯变形，包括以下步骤：
步骤4-1、对第n+1步的六面体单元进行逆向平移：以第n+1步的六面体单元质心为参考点，整体逆向平移至第n步的六面体单元质心上；
步骤4-2、对步骤4-1中经过逆向平移的六面体单元进行第一次逆向转动：将经过逆向平移的六面体单元水平方向上的四个平面质心构造参考平面求解法向量，同理，对第n步的六面体单元也构造出参考平面法向量，以两个法向量的垂直向量为转动向量，夹角为转动角，经过逆向平移的六面体单元质心为转动原点进行整体逆向转动；
步骤4-3、对步骤4-2中经过第一次逆向转动的六面体单元进行第二次逆向转动：将经过第一次逆向转动的六面体单元垂直方向上的四个平面质心构造参考平面求解法向量，同理，对第n步的六面体单元也构造出参考平面法向量，以两个法向量的垂直向量为转动向量，夹角为逆向转动角，经过第一次逆向转动的六面体单元质心为转动原点进行整体逆向转动；
将经过第二次逆向转动的六面体单元质点位置与第n步的六面体单元质点位置作差，计算出质点纯变形：
Figure FDA0002964108420000013
式中，qi为六面体单元各质点的纯变形；
Figure FDA0002964108420000014
为第n步的各质点位置；xi″′为经过逆向运动后的各质点位置；
步骤5、计算质点内力：以第n步六面体单元质心为原点建立局部坐标系
Figure FDA0002964108420000021
将全局坐标系O-xyz下各个质点位置和纯变形转换到局部坐标系中，利用单元形函数求解出虚拟质点内力，将虚拟质点内力正向运动，获得实际位置的质点内力；
步骤6、根据质点中央差分运动公式计算质点位移xn+1；
步骤7、判断分析步数是否达到总步数，当n+1＝N时，达到所需总步数，循环计算结束，否则转达步骤3。
2.根据权利要求1所述的一种六面体单元有限质点法，其特征在于，步骤3更新质点位移的方式为：将第n步和第n-1步的质点位移xn和xn-1作为计算质点位移xn+1的参数。
3.根据权利要求1所述的一种六面体单元有限质点法，其特征在于，步骤6计算质点位移xn+1的中央差分运动公式为：
Figure FDA0002964108420000022
式中，Pn＝fn+Fn，为第n步质点外力和；fn为第n步质点内力；Fn为第n步质点外力。
CN202110246003.3A 2021-03-05 2021-03-05 一种六面体单元有限质点法 Pending CN113032974A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110246003.3A CN113032974A (zh) 2021-03-05 2021-03-05 一种六面体单元有限质点法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110246003.3A CN113032974A (zh) 2021-03-05 2021-03-05 一种六面体单元有限质点法

Publications (1)

Publication Number Publication Date
CN113032974A true CN113032974A (zh) 2021-06-25

Family

ID=76468511

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110246003.3A Pending CN113032974A (zh) 2021-03-05 2021-03-05 一种六面体单元有限质点法

Country Status (1)

Country Link
CN (1) CN113032974A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113849993A (zh) * 2021-08-31 2021-12-28 清华大学 一种基于机械臂柔性变形的工业机器人动力学建模方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2017037553A (ja) * 2015-08-12 2017-02-16 国立研究開発法人情報通信研究機構 運動解析装置および運動解析方法
JP2019175213A (ja) * 2018-03-29 2019-10-10 国立研究開発法人情報通信研究機構 モデル解析装置、モデル解析方法、およびモデル解析プログラム
CN110750933A (zh) * 2019-11-19 2020-02-04 北京理工大学 一种耦合Lagrange质点和Euler方法的精确界面追踪处理方法
WO2020242244A1 (ko) * 2019-05-30 2020-12-03 엘지전자 주식회사 포인트 클라우드 데이터 처리 방법 및 장치

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2017037553A (ja) * 2015-08-12 2017-02-16 国立研究開発法人情報通信研究機構 運動解析装置および運動解析方法
JP2019175213A (ja) * 2018-03-29 2019-10-10 国立研究開発法人情報通信研究機構 モデル解析装置、モデル解析方法、およびモデル解析プログラム
WO2020242244A1 (ko) * 2019-05-30 2020-12-03 엘지전자 주식회사 포인트 클라우드 데이터 처리 방법 및 장치
CN110750933A (zh) * 2019-11-19 2020-02-04 北京理工大学 一种耦合Lagrange质点和Euler方法的精确界面追踪处理方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
李宏光等: "基于六面体体积坐标的新型8结点实体单元", 清华大学学报(自然科学版), vol. 49, no. 11, 15 November 2009 (2009-11-15), pages 1856 - 1860 *
王爱民等: "有限元计算中疏密网格间过渡单元的构造", 清华大学学报(自然科学版), vol. 39, no. 8, 10 August 1999 (1999-08-10), pages 100 - 103 *
王震等: "基于六面体网格的向量式有限元分析及应用", 计算力学学报, vol. 35, no. 4, 31 August 2018 (2018-08-31), pages 480 - 486 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113849993A (zh) * 2021-08-31 2021-12-28 清华大学 一种基于机械臂柔性变形的工业机器人动力学建模方法
CN113849993B (zh) * 2021-08-31 2022-10-25 清华大学 一种基于机械臂柔性变形的工业机器人动力学建模方法

Similar Documents

Publication Publication Date Title
Jaiman et al. Assessment of conservative load transfer for fluid–solid interface with non‐matching meshes
Zhang et al. Analysis of fluid–structure interaction problems with structural buckling and large domain changes by ALE finite element method
CN109726465B (zh) 基于非结构曲边网格的三维无粘低速绕流的数值模拟方法
CN107984472A (zh) 一种用于冗余度机械臂运动规划的变参神经求解器设计方法
CN112016167A (zh) 基于仿真和优化耦合的飞行器气动外形设计方法及系统
CN112001004B (zh) 一种解析中高频振动结构能量密度场的nurbs等几何分析方法
Sadeghi et al. Application of three-dimensional interfaces for data transfer in aeroelastic computations
CN113032974A (zh) 一种六面体单元有限质点法
Fernández Entropy-stable hybridized discontinuous Galerkin methods for large-eddy simulation of transitional and turbulent flows
CN113128085B (zh) 一种基于状态空间形式涡格法的阵风减缓分析方法
CN111341449A (zh) 一种虚拟血管介入手术训练的模拟方法
CN112580240B (zh) 一种适用于复杂大柔性飞机建模的非线性子结构方法
Yu et al. RBFs-MSA hybrid method for mesh deformation
CN113378440A (zh) 一种流固耦合数值模拟计算方法、装置及设备
Strofylas et al. An agglomeration strategy for accelerating RBF-based mesh deformation
CN110705154B (zh) 航空器开环气动伺服弹性系统模型均衡降阶的优选方法
CN117390778A (zh) 一种考虑薄板应变强化效应的塑性安定上限载荷计算方法
Rong et al. Hybrid finite element transfer matrix method and its parallel solution for fast calculation of large-scale structural eigenproblem
CN112949104A (zh) 一种协作机器人实时模态分析方法
Canfield et al. Shape continuum sensitivity analysis using Astros and caps
Bogomolov et al. Comparative verification of numerical methods involving the discontinuous shapeless particle method
CN110674599B (zh) 气动伺服弹性系统非定常气动载荷的有理近似优化方法
Cella et al. Development and validation of numerical tools for FSI analysis and structural optimization: The EU RIBES project status
Xia et al. Analysis-aware modelling of spacial curve for isogeometric analysis of Timoshenko beam
Drofelnik et al. Rapid calculation of unsteady aircraft loads

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