MeshFree-技术原理

midasMeshfree Analysis Manual

摘要：

      传统的有限元分析,需要几何模型的简化和网格的划分,模型越复杂,模型简化工作以及网格划分对工程师的技术和经验要求就越高,这无疑会给广大工程师造成一定的困扰。迈达斯技术有限公司自1989年研发CAE，在实体仿真领域,开发了基于隐式边界法(IBM,Implicit boundary method)的分析软件midas MeshFree，它不需要进行模型的简化和网格的划分，只需导入几何模型、添加荷载和边界条件以及后处理三个分析步骤，降低了工程师的使用难度，可以大大地缩短产品开发周期，提高生产效率。

关键词：CAE；隐式边界法 ( IBM）；midas MeshFree；无网格划分

 D函数的构造   

     midas MeshFree用的是隐式边界 法(IBM)[1,2]。简便起见，以线性静力分析为例进行说明。如图1-1所示的求解问题， 首先在分析对象上生成结构化网格 [3]，其单元非常的规则。

这些单元可以分为三大类，一种是完全在分析对象内部的，图 1-1粉红色线框部分，称为内部单元；第二种是在分析对象边界上的，称为边界单元 。边界单元又可以分两种，一种是含有边界条件的，图 1-1蓝色线框部分，一种是不含边界条件的 ，图1-1绿色线框部分，第三种是完全在分析对象外部的，叫外部单元， 图1-1紫色线框部分，实际上它不需要计算。构造位移函数如下

D是 Dirichlet函数，N是形函数，为网格上分片插值/逼近的网格变量，是节点强制位移。这样构造函数的目的是为了满足边界条件。D函数的表达式如下:

(3)式中 D函数里面的 ϕ函数是用来描述边界的，是隐式方程。当 ϕ≤0时， D函数的值为0，表示在边界上或分析对象外面；当 0≤ϕ≤δ时 δ 为较小的常数一般或者更小,表示在靠近边界的窄带范围内；ϕ≥δ时，表示在分析对象内部。

刚度矩阵的计算

虚功方程：

上式中，{σ}表示应力张量，{t}表示面向量， {b} 表示体力向量，代入 (2)式：

对于内部单元，因为 D函数的值为1,得到的刚度矩阵表达式如下：

(4)式与传统有限元法得到的式子一样.因此,对于内部单元,采用通用的有限元方法处理。

对于不含有边界条件的边界单元,D函数的值也为1。这类单元被分为有材料和无材料两部分，通过调整高斯积分点的权重系数和高斯积分点的位置来实现单元刚度矩阵的计算。因而刚度矩阵的形式与式 (4)一致。

对于含有边界条件的边界单元，位移表达式中含有 D函数，所以应变矩阵中含有D函数的导数。因为δ的值很小，所以D函数及其导数对刚度矩阵的影响不可忽略。此时，刚度矩阵的表达式如下：

对式(5)中的每一项进行积分 ，并相加得到含有边界条件的边界单元刚度矩阵。

计算每个单元的刚度矩阵以及荷载向量，并进行组装，最后求解线性代数方程组得到结果。从上面大致的推导过程可以看出，MeshFree所采用的方法与有限单元法是很类似的，所以它的计算精度是有保障的。

边界值函数的构造

边界值函数的作用是使构造的位移函数满足强制位移，用形函数来构造位移函数，能够保证解结构的完备性要求。

上式中，是单元的形函数，是单元第 i个节点的强制位移。

在结构的边界上施加强制位移，需要在边界经过的所有单元节点上施加强制位移强制位移，其它不相关的节点施加强制位移0。

1. Introduction

   This manual describes theconcepts and theories employed in midas Meshfree. The material presented in themanual is aimed at users with basic knowledge in computational mechanics andrelated engineering fields. As many users are familiar with finite elementmethod (FEM), descriptions are made extensively utilizing notations andalgorithm that originate from classical FEM.

1.1  Introduction to Meshfreeconcept

Although finite element technology andcommercial applications have matured and proven in industrial usage to boostproductivity, small businesses still lag behind in adopting simulationtechnologies in product designs. In this regard, the immediate obstacle isconverting the design CAD model into finite element discretized analysis model.

The “meshless” and “meshfree” terminologies incomputational mechanics have been around with abundant researches conducted.Some methods are targeted at handling extreme nonlinearities while others areadequate for multiscale modeling. The “meshfree” concept in midas Meshfreeemphasizes the ease of use and is directly targeted at product designers; beingfree from laborious process of mesh generation and related cleanup.

1.2 Notations and signconventions

midas Meshfree uses the followingsign conventions for stress for consistency and to avoid confusion. Note thatpositive normal stress refers to tension state while negative normal stress thecompression state.

2.  Basic concepts andformulation

This chapter describes the basic concepts andformulation of midas Meshfree. The key concept that is unique compared totraditional FEM can be summarized as follows.

-The mesh or grids do not conform tothegeometry

The domain of computation is immersed insidestructured grids. Size of grids can be controlled independently for individualCAD part. Similar concepts can be found in finite cell method (ref. 2.1),immersed boundary method in the computational fluid dynamics (ref 2.2).

-Grids are automatically generated bythesolver

For given three-dimensional CAD model, thegrids are generated automatically without user intervention. Unlike automaticmesh generation in FEM which fail more often than not without geometrysimplifications and clean ups, success rate is nearly one hundred percent.

2.1  Compute domain immersed grid elementformulation

For linear elasticity, Hu-Washizu variationalprinciple (ref 2.3, 2.4) can be regarded as the most fundamental equationencompassing equations of equilibrium, constitutive and compatibility:

This generic form can be reduced toHellinger-Reissner variational principle (ref 2.5, 2.6) under the assumptionthat strain and stress satisfy the constitutive equation.

This can further betransformed to the principle of virtual work when compatibility relations between εand Ñ_eu are satisfied.For the sake of simplicity, further derivationassumes linear elastic structural analysis based on assumed displacement,single field formulation.

Figure 2.1shows spatial grid withcomputational domain, Wimmersed inside. The computational domain is discretized by gridelements in a nonconforming manner. Each grid element can be classified intointerior and boundary grid elements. Boundary grid elements can be furtherdivided into those with Dirichlet boundary condition applied and those without,i.e. withNeumann boundary condition. For convenience, interior and Neumann boundary grid elementdomain is designated asand Dirichlet boundary grid element domainis designated as

Similar to finite element procedure,we need to restrict the integral domain to single grid element for applyingeither of the variational principles. For interior grid elements and boundaryelements with Neumann boundary condition, the displacement vector, u is interpolated byshape functions and associated nodal degrees of freedom values.

For boundary grid elements whereDirichlet boundary condition of is prescribed, equation [2-4] is modified in a way that theboundary condition is naturally satisfied at S_u

where  is a diagonal matrix ofspatial functions that modifies solution structure. In essence together withprescribed solution vector, this function enables the solution to satisfy theDirichlet boundary condition at nonconforming surface:

Numerically, the special function is adoptedsuch that the function value changes spatially from one to zero near the vicinity of theDirichlet boundary surface, S_u

Based on equations [2-4] and [2-5], the correspondinggrid element strains can be represented as,Principle of virtual work in equation [2-3] canbe restated in terms of internal and external force vector in the followingform.

Employingstrain-displacement relations expressed equations [2-7] and [2-8], the globalinternal force vector is expressed in terms of summation of contributions fromeach grid element as,where each grid element stiffnessmatrices and force vectors are expressed in the following form.2.2 Other kinematicconstraints

In finite element method, kinematic constraints are usually described in terms of nodal variables and incorporated into the global system matrix as multi-point constraints. However, in midas Meshfree, kinematic constraints are defined over geometries where nodes are not present. Thus, special displacement fields that inherently satisfy the kinematic constraints are devised and incorporated. Among the kinematic constraints available, rigid body and welded contact constraints are introduced herein

Considering two surfaces separated by adistance constrained under welded contact condition, with additional rotationalmotion considered, the geometric nonlinear kinematic relations between two corresponding points in contact surfaces are represented as,

2.4 Numerical integration

   The key to the accuracy of numerical methods based on geometrically non-conforming grids is the domain integration method where the domain is provided as arbitrary threedimensional  CAD data. Robustness, accuracy as well as computational efficiency highly depend on volume and surface integration techniques.

For volume integration, several methods are available in midas Meshfree:

1) Gauss quadrature method with recursive octree domain division.

Figure 2.4 (a) shows the grid element division concept in 2D where the domain is divided along the domain boundary recursively. The recursion continues until integration error is within acceptable tolerance.

2) Gauss quadrature method with domain division by tetrahedrons.

As shown in Figure 2.4 (b), the domain defined by intersection of computational domain and of grid element domain

is subdivided into tetrahedrons. These tetrahedrons are used as base geometry for integration purpose only.

3) Proprietary method based on adaptive generation of integration points aided by raytracing algorithm with automatic geometry repair. The adaptive procedure is driven by minimizing integration error as well as reducing the overall number of integration points.

By default, volume integration is conducted using the third method. It can repair the flaws in CAD models which in many cases exist. Other methods are utilized to make the whole integration procedure fail-safe.

With regards to surface integration method, the CAD surface is divided into triangular patches, and integration is conducted using Gauss quadrature on triangles

2.5 Hybridformulation

Finite elements especially low-order types based on assumed displacement formulation generally suffer from poor accuracy due to locking phenomenon. This fact extends to the immersed grid method adopted in midas Meshfree. To enhance accuracy mixed hybrid formulation is developed based on Hellinger-Reissner variational principle (see equation [2-2]).

In addition to the displacement assumption in equations [2-4] and [2.5], stress is independently assumed in the grid element level (ref 2.7):For interior grid elements and boundaryelements with Neumann boundary condition, the right hand side of equation [2-2]can be expressed as,For boundary grid element with Dirichletboundary condition prescribed, the same terms takes the following form:

2.6 Error estimation

Error is introduced during the discretizationprocess of the solution field into finite number of grid elements. To assessthe discretization error, a posteriori errorestimation can be utilized. To this end, strain energy error is calculatedbased on a stress error definition. The stress error can represented as errorbetween exact stress that satisfies the governing equation and approximate gridstress which includes errors fromdiscretization.Similar to the methods used in the finiteelement method community, stress obtained by a separate recovery procedure isemployed as an estimate to the exact stress to define stress error.Meshfree provides error measures as well asstrain energy for individual part in an assembly, to assess discretizationerror. The information can be useful for adjusting overall grid size or partlevel grid size. Table 2.1 summarizes the output information generated toassess discretization error.

Table 2.1 Output information to assess discretization error

2.7  Heat transfer analysis using immersed grid elements

The heat transfer formulation is based on a variational form of the energy balance equation with Fourier law of heat conduction.Loading and boundary conditions includeprescribed temperature, heat generation on volume, surface flux, surfaceconvection and radiation.For transient heat transfer analysis, backwarddifference algorithm is adopted. The final integral equation is expressed as:

If material property is dependent ontemperature or in the presence of radiation, the equation is nonlinear. Thenonlinear equation is solved using Newton-Raphson method, where linearizedresidual equation is iteratively solved until convergence.

2.7 Reference

1. 2.1 A. Duster, J Parvizian, Z Yang, andE Rank, “The finite cell method for three-dimensional problems of solidmechanics”, The finite cell method for three-dimensional problems of solidmechanics, Vo. 197, Issues 45-48,2008

2. 2.2 R Mittal and G Iaccarino, “Immersedboundary methods”, Annual Review of Fluid Mechanics, Vo. 37, 2005

3. 2.3 Hu,H.C., “On some variational principles in the theory of elasticity and thetheory of plasticity,” Scintia Sinica, Vol. 4,1955

4. 2.4 Washizu,K., On the Variational Principles of Elasticity, Aeroelastic and StructuralResearch Laboratory, MIT, Technical Report,1955

5. 2.5 Hellinger,E., “Der allgemeine Ansatz der Mechanik der Kontinua,” Encyclopadie derMathernafischen Wissenschaften, Vol. 4,1914

6. 2.6 Reissner,E., “On a variational theorem in elasticity,” Journal of Mathematical Physics,Vol. 29, 1950

     2.7 Pian, T.H.H. and Sumihara, K., “Rational approach forassumed stress finite elements,” International Journal for Numerical Methods inEngineering, Vol. 20,1984