CN113343523B - 一种基于空间网格的有限元数值模拟分析方法 - Google Patents

Abstract

本发明涉及一种基于空间网格的有限元数值模拟分析方法，其特征在于，包括以下步骤：将物理场分析计算的空间进行网格分割，分割得到空间网格；将待分析的物体离散化为颗粒；待分析的物体发生形变，则颗粒进行变化，但空间网格不发生变化；对发生形变后的待分析的物体进行下一步分析计算时，空间网格被包含的颗粒重新赋予新的物理量场，对新的物理量场再进行有限元分析计算，得到各个节点的物理值，进而根据各个节点的物理值内插计算得到颗粒的物理值。本发明的物理场分析计算是运用空间网格进行有限元分析；被分析物体的物理值计算，由离散化每一颗颗粒由颗粒所在的空间网格的节点内插计算得到，被分析物体按照得到物理值进行变化。

Description

一种基于空间网格的有限元数值模拟分析方法

技术领域

本发明涉及一种运用有限元的数值模拟分析方法。

背景技术

有限元法，是将分析物体划分成有限的网格，进行物理场分析计算的数值模拟方法。从21世纪六十年代开始，作为一种数值模拟的有效方法，从土木工程到航空飞行，有限元法都得到了广泛地运用。同时，有限元法还运用到金属成型领域中，分析计算金属在变形工程中的应力和应变。但是，传统的有限元法有一个致命的缺点，在分析金属的大变形时，由于网格的不规则大变形，计算经常被迫中断，无法进行下去。

为解决网格的不规则大变形引起的分析计算中断，现有的主要方法是对不规则区域进行网格重新划分，生成新的规则的网格，使得计算能够重新继续。这种方法存在以下问题：

1)失真：计算中途对网格进行了重新划分，失去了网格原始的形状，失去了真实的信息。

2)费时：网格的重新划分花费时间，有时还需要多次的人工干预。

3)局限性：在大变形的分析中，对于复杂的形状，要得到一个新的且可以适用于有限元计算的网格非常困难，无法使计算重新继续。

发明内容

本发明要解决的技术问题是：传统的有限元法在分析金属的大变形时，由于网格的不规则大变形，计算经常被迫中断，无法进行下去。

为了解决上述技术问题，本发明的技术方案是提供了一种基于空间网格的有限元数值模拟分析方法，其特征在于，包括以下步骤：

步骤1、将物理场分析计算的空间进行网格分割，分割得到m个网格，将每个网格定义为空间网格，且每个空间网格具有n个节点，n≥3；

步骤2、将待分析的物体离散化为一个一个的点，将这种点定义为颗粒，颗粒由横竖的线连接，放置在空间网格中；每一颗颗粒的物理值由颗粒所在的空间网格的节点的物理值内插计算得到；

步骤3、待分析的物体发生形变时，则颗粒按照通过所在网格的节点内插得到物理值进行变化，但是，空间网格不发生变化；

步骤4、对发生形变后的待分析的物体进行下一步分析计算时，空间网格被包含的颗粒重新赋予新的物理量场，对新的物理量场再进行有限元分析计算，得到各个节点的物理值，进而根据各个节点的物理值内插计算得到颗粒的物理值。

此方法用于对在模具作用下发生形变的物体进行有限元数值模拟分析。

，所述物理值为速度，则步骤1中，利用有限元法的公式，对于每个空间网格通过对式(1)所示的变形能量函数Φ最小化得到空间网格的节点的速度：

式(1)中，

表示相当应力，

表示相当应变速度，|τ_f|表示剪切摩擦应力，Δv表示相对速，λ是拉哥朗日系数，V表示体积，S表示面积，

表示空间网格的静水应变速度；

空间网格的静水应变速度

由式(2)定义：

式(2)中，

表示X方向的应变速度，

表示Y方向的应变速度，

表示Z方向的应变速度；

式(1)对于空间网格的节点速度的偏微分公式表示为下式(3)：

式(3)中，ν_i表示节点i处的节点速度，n表示每个空间网格具有的节点的总数，τ_j表示要素j的拉哥朗日系数，m表示空间网格的总数。

优选地，空间网格有4个节点，则每一颗颗粒的速度采用以下步骤计算得到：

步骤a、设颗粒P(x,y)在某个空间网格中，动态速度场的边界条件满足下式(4)：

ul₂+vl₁＝a (4)

式(4)中，l₁、l₂是颗粒P(x,y)在与模具接触点的表面法向矢量的余弦的值，a是模具在此接触点的表面法向速度，u是颗粒P(x,y)沿X轴方向的速度分量，v是颗粒P(x,y)沿Y轴方向的速度分量。

u、v表示为下式(5)：

式(5)中，u₁、u₂、u₃、u₄为四个节点沿X轴方向的速度分量，v₁、v₂、v₃、v₄为四个节点沿Y轴方向的速度分量，N₁、N₂、N₃、N₄为空间网格的四个节点在颗粒处的内插函数的值；

将式(5)代入式(4)，则有：

(N₁u₁+N₂u₂+N₃u₃+N₄u₄)l₂+(N₁v₁+N₂v₂+N₃v₃+N₄v₄)l₁＝a (6)

步骤b、计算颗粒P(x,y)沿X轴方向的速度u以及颗粒P(x,y)沿Y轴方向的速度v，则有：

u＝a₀+a₁x+a₂y+a₃xy＝{1 x y xy}{a}＝{1 x y xy}[A]^-1{u} (7)

v＝b₀+b₁x+b₂y+b₃xy＝{1 x y xy}{b}＝{1 x y xy}[A]^-1{v} (8)

式(7)及式(8)中，

系数a₀、a₁、a₂、a₃以及系数b₀、b₁、b₂、b₃是和节点坐标值有关的值；

(x₁,y₁)、(x₂,y₂)、(x₃,y₃)、(x₄,y₄)是四个节点的坐标，[A]^-1是[A]的逆矩阵；

(x,y)是颗粒P(x,y)的坐标值。

与现有技术相比，本发明具有如下特点：

1)空间网格的运用

物理场分析计算是运用空间网格进行有限元分析；被分析物体的物理值计算，由离散化每一颗颗粒由颗粒所在的空间网格的节点内插计算得到，被分析物体按照得到物理值进行变化。

2)空间网格的边界条件的设定方法。

本方案中，物体(颗粒)的边界条件(比如：速度、温度等)通过内插函数赋加到空间网格的节点上，使得边界条件更准确，更方便地反映到有限元计算分析中，同时，减少了计算对空间网格形状，位置的要求，使用更加方便。

3)空间网格的便捷

空间网格，减少了有限元时操作的复杂程度，使得更多不同教育背景的人都可以理解使用有限元。

4)空间网格的广泛性

空间网格可以方便地运用到各种有限元分析，比如：弹性变形分析、温度分布分析等线性和非线性物理场的计算分析。

附图说明

图1a为模具和物体(颗粒)示意图；

图1b为空间网格示意图；

图1c为空间网格和物体(颗粒)的关系示意图；

图2为边界条件的设置示意图；

图3为应用空间网格计算的大变形结果；

图4为物体(颗粒)的各压缩量时的变化示意图，其中，(a)、(b)、(c)、(d)、(f)分别表示压缩量分别为0、15.4、20.2、21.6、24.0时的变化示意图。

具体实施方式

下面结合具体实施例，进一步阐述本发明。应理解，这些实施例仅用于说明本发明而不用于限制本发明的范围。此外应理解，在阅读了本发明讲授的内容之后，本领域技术人员可以对本发明作各种改动或修改，这些等价形式同样落于本发明所附权利要求书所限定的范围。

本发明提出一种基于空间网格法的有限元法，如图1a、图1b及图1c所示，在本发明提供的方法中，将物理场分析计算的空间进行网格分割，这种网格被定义为“空间网格”，每个空间网格一般具有四个节点。物理场分析计算是运用空间网格进行有限元分析，将待分析的物体进行离散化，分成一个一个的点，将这种点定义为“颗粒”。颗粒由横竖的线连接，放置在网格空间中。每一颗颗粒的物理值由颗粒所在的空间网格的节点内插计算得到。待分析的物体在模具作用下发生形变后，颗粒按照得到物理值进行变化，但是，空间网格不发生变化。待分析的物体发生形变后，进行下一步分析计算时，空间网格被包含的颗粒重新赋予新的物理量场，对新的物理量场再进行有限元分析计算，得到各个节点的物理值。根据内插计算得到颗粒的物理值，进行新的一步变化。

本发明中，空间网格用于节点的物理值计算，再通过内插计算得到颗粒的物理值，由颗粒的物理值表示物体的变化。所以，在计算过程中，运用有限空元计算的空间网格的形状和位置不发生变化，保持原始规则的形状，不会引起计算的中断，使得金属成形的大变形分析计算能够实现。

在上述步骤中，物理值为速度，则利用有限元法的公式，对于每个空间网格通过对式(1)所示的变形能量函数Φ最小化可以得到空间网格的节点的速度：

式(1)中，

表示相当应力，

表示相当应变速度，|τ_f|表示剪切摩擦应力，Δv表示相对速，λ是拉哥朗日系数，V表示体积，S面积表示，

表示空间网格的静水应变速度。

空间网格的静水应变速度

由式(2)定义：

式(2)中，

表示X方向应变速度，

表示Y方向应变速度，

表示Z方向应变速度表示。

式(1)对于空间网格的节点速度的偏微分公式表示为下式(3)：

式(3)中，ν_i表示节点i处的速度(，n表示每个空间网格具有的节点的总数，τ_j表示要素j的拉哥朗日系数，m表示空间网格的总数。

基于计算得到的空间网格的节点的速度采用以下步骤计算得到颗粒的速度：

步骤1、如图2所示，颗粒P(x,y)在某个空间网格中，序号1至4表示该空间网格的4个节点，4个节点和模具表面接触。为了保证颗粒P(x,y)沿模具表面移动，动态速度场的边界条件必须满足下式(4)：

ul₂+vl₁＝a (4)

式(4)中，l₁、l₂是颗粒在与模具接触点的表面法向矢量的余弦的值，a是模具在此接触点的表面法向速度，u是颗粒沿X轴方向的速度分量，v是颗粒沿Y轴方向的速度分量。

u、v表示为下式(5)：

式(5)中，u₁、u₂、u₃、u₄为通过上文计算得到的如图2所示的四个节点沿X轴方向的速度分量，v₁、v₂、v₃、v₄为通过上文计算得到的如图2所示的四个节点沿Y轴方向的速度分量，N₁、N₂、N₃、N₄为空间网格在颗粒处的内插函数的值。

将式(5)代入式(4)，则有：

(N₁u₁+N₂u₂+N₃u₃+N₄u₄)l₂+(N₁v₁+N₂v₂+N₃v₃+N₄v₄)l₁＝a (6)

步骤2、计算利用模具加工的物体的变形速度，包括以下内容：

加工物体的变形速度可以通过所在空间网格的节点的速度计算得到：

式(7)中，系数a₀、a₁、a₂、a₃以及系数b₀、b₁、b₂、b₃是和颗粒P(x,y)的坐标值有关的值，由以下方法确定：

式(8)及式(9)中，(x₁,y₁)、(x₂,y₂)、(x₃,y₃)、(x₄,y₄)是如图2所示的4个节点的坐标。

则有：

[A]{a}＝{u} (10)

[A]{b}＝{v} (11)

式(10)及式(11)中，

则有：

{a}＝[A]^-1{u} (12)

{b}＝[A]^-1{v} (13)

式(10)及式(11)中，[A]^-1是[A]的逆矩阵，是空间网格的4个节点的坐标值的函数，则：

加工物体的某个颗粒沿X轴方向的速度u表示为：

u＝a₀+a₁x+a₂y+a₃xy＝{1 x y xy}{a}＝{1 x y xy}[A]^-1{u} (14)

加工物体的某个颗粒沿Y轴方向的速度v表示为：

v＝b₀+b₁x+b₂y+b₃xy＝{1 x y xy}{b}＝{1 x y xy}[A]^-1{v} (15)

因此，加工物体颗粒的应变速度可以用式(14)和式(15)计算。例如：对平面应变的加工变形，应变速度可以计算如下：

则应力{σ}用下式计算：

式(17)中，[D]是刚性矩阵。

在如图3所示的计算中，物体在高度方向被压缩了57％的情况下，本方法仍然可以有效地进行大变形的计算和分析，而不需要进行网格的重新分割。

如图4所示，在发生大变形时，由物体颗粒组成的4边形形状，被压缩成一个扁平的形状，传统的有限元法的计算是无法继续进行下去。本方法运用了空间网格进行有限元计算和分析，不受到物体大变形的影响而发生计算中断。

Claims (2)

1.一种基于空间网格的有限元数值模拟分析方法，其特征在于，用于对在模具作用下发生形变的物体进行有限元数值模拟分析，设物理值为速度，所述有限元数值模拟分析方法包括以下步骤：

步骤1、将物理场分析计算的空间进行网格分割，分割得到m个网格，将每个网格定义为空间网格，且每个空间网格具有n个节点，n≥3，包括以下步骤：

利用有限元法的公式，对于每个空间网格通过对式(1)所示的变形能量函数Φ最小化得到空间网格的节点的速度：

式(1)中，

表示相当应力，

表示相当应变速度，|τ_f|表示剪切摩擦应力，Δv表示相对速，λ是拉哥朗日系数，V表示体积，S表示面积，

表示空间网格的静水应变速度；

空间网格的静水应变速度

由式(2)定义：

式(2)中，

表示X方向应变速度，

表示Y方向应变速度

表示Z方向应变速度；

式(1)对于空间网格的节点速度的偏微分公式表示为下式(3)：

式(3)中，ν_i表示节点i处的速度，n表示每个空间网格具有的节点的总数，τ_j表示要素j的拉哥朗日系数，m表示空间网格的总数；

步骤2、将待分析的物体离散化为一个一个的点，将这种点定义为颗粒，颗粒由横竖的线连接，放置在空间网格中；每一颗颗粒的物理值由颗粒所在的空间网格的节点的物理值内插计算得到；

步骤3、待分析的物体发生形变，则颗粒按照得到物理值进行变化，但是，空间网格不发生变化；

步骤4、对发生形变后的待分析的物体进行下一步分析计算时，空间网格被包含的颗粒重新赋予新的物理量场，对新的物理量场再进行有限元分析计算，得到各个节点的物理值，进而根据各个节点的物理值内插计算得到颗粒的物理值。

2.如权利要求1所述的一种基于空间网格的有限元数值模拟分析方法，其特征在于，空间网格有4个节点，则每一颗颗粒的速度采用以下步骤计算得到：

步骤a、设颗粒P(x,y)在某个空间网格中，动态速度场的边界条件满足下式(4)：

ul₂+vl₁＝a (4)

式(4)中，l₁、l₂是颗粒P(x,y)在与模具接触点的表面法向矢量的余弦的值，a是模具在此接触点的表面法向速度，u是颗粒P(x,y)沿X轴方向的速度分量，v是颗粒P(x,y)沿Y轴方向的速度分量；

u、v表示为下式(5)：

式(5)中，u₁、u₂、u₃、u₄为四个节点沿X轴方向的速度分量，v₁、v₂、v₃、v₄为四个节点沿Y轴方向的速度分量，N₁、N₂、N₃、N₄为空间网格的四个节点在颗粒处的内插函数的值；

将式(5)代入式(4)，则有：

(N₁u₁+N₂u₂+N₃u₃+N₄u₄)l₂+(N₁v₁+N₂v₂+N₃v₃+N₄v₄)l₁＝a (6)

步骤b、计算颗粒P(x,y)沿X轴方向的速度u以及颗粒P(x,y)沿Y轴方向的速度v，则有：

u＝a₀+a₁x+a₂y+a₃xy＝{1 x y xy}{a}＝{1 x y xy}[A]^-1{u} (7)

v＝b₀+b₁x+b₂y+b₃xy＝{1 x y xy}{b}＝{1 x y xy}[A]^-1{v} (8)

式(7)及式(8)中，

系数a₀、a₁、a₂、a₃以及系数b₀、b₁、b₂、b₃是和节点坐标值有关的值；

(x₁,y₁)、(x₂,y₂)、(x₃,y₃)、(x₄,y₄)是四个节点的坐标，[A]^-1是[A]的逆矩阵；

(x,y)是颗粒P(x,y)的坐标值。

Images