高斯积分点在有限元中应用

高斯公式应用案例高斯公式是数学中一个非常重要的公式，它在很多领域都有广泛的应用。

本文将介绍几个关于高斯公式应用的案例，分别来自物理学、工程学和金融学领域。

物理学：电场中的高斯定律高斯公式最早是由德国数学家高斯提出的，但在物理学中也有广泛的应用。

电场中的高斯定律就是一个非常经典的例子。

根据高斯定律，通过一个闭合曲面的电通量等于该闭合曲面内的电荷总量的1/ε0倍，其中ε0为真空介电常数。

这个定律在物理学中被广泛用于计算电场的分布。

我们可以通过高斯定律来计算一个均匀带电球的电场分布，或者通过选择适当的高斯曲面来计算复杂形状电荷分布的电场。

通过高斯公式的应用，我们可以更好地理解电场的性质，对电磁学的学习和实践有很大的帮助。

工程学：有限元分析中的面积分在工程学中，高斯公式也有着举足轻重的地位。

有限元分析是工程学领域中常用的一种数值分析方法，用于求解复杂结构的应力、位移和变形等问题。

在有限元分析中，经常需要对复杂的形状进行面积分计算，而高斯公式可以帮助我们高效地进行这些积分。

通过高斯公式，我们可以将复杂形状的面积分转化为一系列关于标准形状的积分，从而更方便地进行数值计算。

这种方法既可以提高计算效率，也可以提高计算的精度，因此在工程学中有着广泛的应用。

金融学：期权定价中的黑-斯科尔斯模型除了自然科学和工程学领域，高斯公式在金融学中也有一些重要的应用。

其中一个著名的例子就是在期权定价中的黑-斯科尔斯模型。

黑-斯科尔斯模型是用于计算欧式期权价格的数学模型，它可以根据标的资产的价格波动情况、执行价格、无风险利率等因素来估算期权的价格。

在这个模型中，高斯公式被用来计算标的资产价格的概率分布。

通过高斯公式，我们可以更准确地估算出期权的价格，对投资者和金融机构来说都具有重要的意义。

通过以上三个领域的案例，我们可以看到高斯公式在自然科学、工程学和金融学中都有着广泛的应用。

它不仅是一个重要的数学工具，也是连接数学与实际应用的桥梁。

各单元高斯点在有限元分析中，高斯点是一种用于数值计算的重要概念。

在分析结构或材料时，通常会将其离散成有限个单元来进行计算。

而每个单元都会包含若干个高斯点，用来近似表示单元内的物理量。

在这篇文章中，我将讨论各单元高斯点的重要性及其在有限元分析中的应用。

首先，让我们来了解一下高斯点的概念。

高斯点是一种数值积分点，通过在这些点上对被积函数进行数值积分，可以得到近似的积分值。

在有限元分析中，高斯点通常被用来近似表示单元内的应力、应变、位移等物理量。

通过在高斯点上对这些物理量进行插值，可以得到整个单元内的近似解。

在有限元分析中，单元是构成整个结构或材料的基本单位。

不同类型的单元可以用来表示不同的几何形状，比如一维线元、二维三角形元、四边形元等。

而每个单元内都包含若干个高斯点，这些高斯点的数量通常是提前确定的。

在进行有限元计算时，我们需要在每个单元内对物理量进行数值积分，而高斯点就是用来指定积分的位置和权重的。

在有限元分析中，各单元高斯点的选择对计算结果的精度和稳定性有着重要的影响。

一般来说，高斯点的数量越多，计算结果越精确，但计算量也会增加。

因此，在实际应用中，需要在计算精度和计算效率之间进行权衡。

通常情况下，工程师会根据具体问题的需求和计算资源的限制来选择合适的高斯点数量。

除了高斯点的数量，高斯点的位置也是至关重要的。

在有限元分析中，高斯点的位置通常被选取在单元的几何中点或者重心附近。

这样可以保证在近似表示单元内物理量时，能够更好地反映单元内部的特性。

而如果高斯点的位置选择不当，可能会导致计算结果出现偏差，甚至无法收敛。

在实际工程中，有限元分析是一种非常重要的工具，可以用来分析结构的应力、变形等物理量。

而各单元高斯点的选择是有限元分析中一个至关重要的环节。

通过合理选择高斯点的数量和位置，可以提高计算结果的精度和可靠性。

因此，在进行有限元分析时，工程师需要仔细考虑各单元高斯点的选择，以确保得到准确的计算结果。

c3d10的积分点C3D10元素是有限元分析软件中常用的一个十节点立方体元素（brick element）。

它对于复杂的三维结构的有限元分析提供了很大的便利，能够有效地模拟力学行为并预测结构的响应。

本文将详细介绍C3D10元素的积分点，以及如何使用这些积分点进行分析。

首先，我们来了解C3D10元素的基本信息。

C3D10元素有十个节点，它按照正立方体的八个顶点和中心点的顺序进行编号。

同时，它有六个面，每个面都有四个节点。

此外，C3D10元素的两条棱上各带有一个节点，总共有十条棱。

这种节点的分布方式使得C3D10元素能够更好地适应结构的几何形状。

在有限元分析中，C3D10元素经常用于模拟具有复杂几何形状的结构，例如汽车车身、飞机机身和建筑物等。

其三维形状可以更好地模拟真实结构，并且由于具有足够的节点数目，能够精确地捕捉结构的力学行为。

在C3D10元素中，为了进行数值积分和力学计算，需要在元素内部选择一系列积分点。

这些积分点通常按照高斯积分点的规则进行选择，以获得更准确的结果。

积分点的数量和位置可以根据具体的分析要求进行选择。

一般来说，积分点越多，结果越准确，但计算量也会增加。

C3D10元素的常见积分点数目有1、8、27和64等。

1个积分点用于代表整个单元，它的位置位于单元的质心上。

8个积分点位于单元的八个顶点上，27个积分点则按照某种规律分布在元素内部的各个位置，64个积分点则进一步增加了积分精度，用于需要更高精度的分析场合。

使用C3D10元素进行有限元分析时，我们需要在每个元素内部选择合适数量的积分点，并计算每个积分点处的场量值。

这些场量通常包括位移、应力、应变等。

通过选择适当的积分点数量，我们能够获得足够准确的场量值，从而在整个结构中分析出力学响应。

C3D10元素的积分点选择与力学计算有着紧密的联系。

在弹性力学中，通常使用高斯积分点进行数值计算，以获得足够精确的结果。

在非线性力学中，由于材料性质的变化，可能需要更多的积分点来获得准确的结果。

二阶六面体单元积分点个数1. 介绍在有限元分析中，常常需要对六面体单元进行积分以计算其质量矩阵、刚度矩阵等。

而对于二阶六面体单元，其积分点的个数是一个重要的参数，它决定了数值计算的精度和效率。

本文将介绍二阶六面体单元积分点个数的计算方法和应用。

2. 二阶六面体单元二阶六面体单元是一种常用的有限元单元，它由六个面构成，每个面都是一个四边形。

它的节点个数为27个，其中8个节点位于六面体的八个顶点上，剩下的19个节点位于六个面的中点上。

二阶六面体单元的形状函数通常采用三次Lagrange插值函数来表示。

3. 积分点个数的计算方法计算二阶六面体单元的积分点个数需要考虑到积分的精度和计算效率。

一般来说，积分点个数越多，计算结果越精确，但计算量也越大。

在实际应用中，需要根据具体问题的要求来选择合适的积分点个数。

一种常用的计算积分点个数的方法是根据六面体单元的阶数和积分点的阶数来确定。

二阶六面体单元的阶数为2，积分点的阶数为3。

根据经验公式，二阶六面体单元的积分点个数可以通过以下公式计算：N = (n+1)^3其中，N为积分点个数，n为积分点的阶数。

对于二阶六面体单元，带入n=3，可以得到N=64，即需要64个积分点来进行积分计算。

4. 积分点个数的应用二阶六面体单元的积分点个数在有限元分析中起着重要的作用。

它决定了计算结果的精度和计算效率。

4.1 精度积分点个数越多，计算结果越精确。

在求解复杂问题时，需要使用更多的积分点来提高计算精度。

通过增加积分点个数，可以减小数值计算误差，得到更接近真实解的结果。

4.2 计算效率积分点个数越多，计算量越大，计算效率越低。

在实际应用中，需要在精度和效率之间进行权衡。

如果问题的精度要求不高，可以选择较少的积分点个数来提高计算效率。

而对于精度要求较高的问题，需要使用更多的积分点来保证计算结果的准确性。

4.3 积分公式二阶六面体单元的积分通常采用高斯积分公式来计算。

高斯积分公式是一种常用的数值积分方法，它通过在积分区间上选择一些特定的节点和权重来近似计算积分结果。

高斯公式应用案例高斯公式是数学上非常重要且广泛应用的公式。

它可以帮助我们计算各种形状的定积分，例如曲线下面积、曲线围成的曲边梯形面积、曲线周长等。

在不同的领域中，高斯公式都有着重要的应用。

本文将介绍高斯公式在物理、工程、经济和生物等领域中的应用案例，以及其在实际问题中的重要性。

一、物理高斯公式在物理学中有着广泛的应用，特别是在电磁学和力学领域。

在电场和磁场中，高斯公式可以用来计算电场线和磁场线的通量，从而求解电荷和磁荷的分布情况。

在引入高斯公式后，可以简化问题求解的复杂度，从而更方便地研究电磁场的性质。

在静电学中，高斯定律描述了电场的产生和分布。

利用高斯公式，可以求解由不同电荷分布所产生的电场强度，进而解决电场环绕导体的分布问题。

高斯公式还可以对电场在不同介质中的分布情况进行精确的描述，为电场的应用提供了重要的理论基础。

在力学领域，高斯公式也常用于计算曲线轨迹下物体的运动状态。

当一个物体沿着曲线运动时，我们可以利用高斯公式来计算曲线上的力的合成，求解物体的加速度和速度等动力学问题。

高斯公式的应用使得复杂的动力学问题变得更加清晰和可计算，为物理学研究提供了重要的数学工具。

二、工程在工程领域，高斯公式也有着重要的应用价值。

特别是在结构分析和流体力学中，高斯公式可以帮助工程师解决各种复杂的结构计算和流体运动问题。

在结构分析中，高斯公式可以用来计算不同形状结构的受力情况。

利用高斯公式可以求解曲线形状的梁在受力作用下的变形和内应力分布，为结构设计提供了重要的数学工具。

在有限元分析中，高斯公式也可以用来建立与结构形状相关的数学模型，进而对结构进行精确的应力分析和应变计算。

在流体力学中，高斯公式被广泛应用于计算流体在不同形状容器中的流动情况。

在管道工程中，高斯公式可以用来计算管道中流体的流速分布和流量情况，从而指导管道的设计和运行。

高斯公式还可以对复杂的流场进行数值模拟和计算，为工程师提供了重要的工具来研究流体动力学问题。

单元节点和积分点有什么区别学过数值积分的应该知道,有限元中的积分点指高斯积分点,因为这些点的收敛性好,精度高。

1、节点7_J1D o(^在单元内,采用形函数来表述单元内变量的分布规律。

而节点值是在节点处的对应物理量。

以简单矩形单元的温度为例:1n-]*}1p"g I }四个节点i,j,m,n的温度分别为Ti,Tj,Tm,Tn.则以单元内自然坐标(x,y),(-1,-1),(-1,1),(1,-1),(1,1)分别为四个节点,单元内温度分布为:T={Si, Sj, Sm, Sn} {Ti, Tj, Tm, Tn}Si=1/4(1-x)(1-y)Sj=1/4(1+x)(1-y)] q v5uSm=1/4(1+x)(1+y) E g@0w4]y6l.?-aSn=1/4(1-x)(1+y)M r S1T~6]+d:C(单元的形函数我们可以从手册中查到)从而我们知道了温度在单元内的分布。

2、积分节点5G z)\ \我们需要对温度在单元内的面积上进行积分时,因为节点的温度显然与x,y无关,我们只需要考虑对形函数积分。

f"X] r1O采用Gauss_Legendre多项式计算积分时,我们只需要计算根据特定积分点的值(在自然坐标系下是固定的,可以查手册,这些点也叫高斯点、积分点)并加以权重就可以。

这就把复杂的积分问题变成了简单的代数问题。

因为形函数只有单元有关,所以积分点也只与单元形状有关。

3.应力一般采用多个积分点的相互插值或外延来计算节点应力。

这只是为了减少误差。

因为在积分点应力比节点具有更高阶的误差。

从理论上说,形函数已知后,用Maple或者Mathematic等软件进行符号积分的话,是可以精确计算出刚度矩阵和质量矩阵,但是这样做的话,对于工程实际应用来说并不合适(9W F;c#j原因:1,费时;2,Mindlin中厚板有剪力锁死问题,有时候需要采用缩聚积分),所以有些书上会把2节点梁单元的刚度阵直接写出来,但是再复杂点的单元,就使用数值积分(Newton-Cotes积分和高斯积分)高斯积分的话,积分点不在节点上9B N V4L2K*o牛顿-科斯的积分点就是节点,这样得到的质量矩阵是集中质量阵形式个人理解:1.节点作用构造形函数,节点的多少描述规则形状单元内的应力的近似分布情况,并获取节点上的位移值m0D8Y.p5p2.积分点作用是构造规则形状单元与曲边(曲面)单元的转化的变换函数,积分点的选取多少和选取的位置直接关系到这种“映射”-y-j A.|K7r的精确程度,刚度矩阵、边界条件的转化都用到了坐标变换的积分关系,一般取高斯积分点能使被积函数计算精度尽量高。

第39卷第3期Vol.39,No.32009年3月J OURNAL OF UNIVE RSITY OF SCIE NCE AND TECHNOLO GY OF CHINAMar.2009文章编号:025322778(2009)0320331206收稿日期:2006211213;修回日期:2007206228基金项目:国家自然科学基金(10632080,10572134)资助.作者简介:王峰,男,1976年生,博士.研究方向:计算力学.E 2mail :wangfeng0706@ 通讯作者:王肖钧,教授.E 2mail :xjwang @冲击力学有限元计算中的数值积分王　峰1,2,王肖钧1,卞　梁1,李建光1,劳　俊1(1.中国科学技术大学近代力学系,安徽合肥230027;2.解放军炮兵学院无人机和机械工程系,安徽合肥230031)摘要:在Lagrange 有限元基础上,介绍了计算运动方程中节点力的不同积分方法.单点积分方法具有较高的计算效率,为了控制沙漏变形,必须引入抗沙漏节点力;采用2×2×2高斯积分可以避免沙漏变形,并有较高的计算精度,但导致计算量增大;而采用局部2×2×2高斯积分则同时具有两者的优点.三维侵彻的计算结果表明局部2×2×2高斯积分能够很好地控制沙漏变形,并有较高的计算效率;一维应变波的模拟计算结果也表明,2×2×2高斯积分比单点积分更加接近理论值.这说明所述方法和所建程序的合理性和有效性,它为侵彻贯穿过程的数值分析提供了一种实用和有效的手段.关键词:爆炸力学;数值模拟;有限元;单点积分;高斯积分中图分类号:O385 文献标识码:ANumerical integration method of f inite elementcomputation in impact mechanicsWAN G Feng 1,2,WAN G Xiao 2jun 1,B IAN Liang 1,L I Jian 2guang 1,L AO J un 1(1.Depart ment of Modern Mechanics ,Universit y of S cience and Technology of China ,Hef ei 230027,China;2.Depart ment of Unmanned A erial V ehicle and Mechanical Engineering ,PL A A rtillery A cadem y ,Hef ei 230031,China )Abstract :Based on t he analysis of Lagrange finite element met hod ,numerical integration met hods to comp ute nodal force were briefly described.In order to cont rol hourglass deformations ,anti 2hourglass nodal forces have to be used for single point quadrat ure ,which has higher efficiency.The hourglass modes would be cont rolled and t he numerical accuracy raised effectively if 2×2×2Gaussian quadrat ure were adapted ,but t his would entail heavy comp utation.However ,partial 2×2×2Gaussian quadrat ure has t he merit s of two kinds of integration met hod.The numerical example of penet ration indicated t hat partial 2×2×2Gaussian quadrat ure could cont rol hourglass deformations effectively and has higher efficiency.And t he numerical simulations of one 2dimensional strain wave also shown t hat 2×2×2Gaussian quadrat ure is much better t han single point quadrat ure.It is concluded t hat t he met hod discussed and t he program we developed are reasonable and effective ,providing a usef ul met hod for t he numerical st udy of penet ration and perforation.K ey w ords :mechanics of explosion ;numerical simulation ;finite element ;single point quadrat ure ;Gaussian quadrat ure0　引言侵彻问题的理论与实验研究始终是新武器研制和高科技领域中的重要课题,无论是进攻或者防御,可靠预测侵彻效果、提高侵彻能力始终是新武器研究中的核心.传统的解析方法以及基于实验研究的经验或半经验方法已经远远不能满足现代高技术武器发展的需要.数值模拟方法,由于具有假设少、物理模型逼真、计算精度高等明显的优势,成为了高速侵彻理论研究中一种主要和广泛认同的方法[1,2].在有限元计算中,采用集中质量法后,守恒方程中的运动方程被离散为通过节点力表示的质量运动方程,节点力的计算需要通过一定的积分方法来实现.为了提高数值模拟侵彻过程的计算效率,通常采用单点积分法,但这种方法的直接后果是立方元被认为是常应力常应变单元,在某些情况下会导致网格相交,甚至网格翻转,计算中断,因此必须通过引入抗沙漏节点力的方法以控制沙漏变形.采用多点积分,如2×2×2高斯积分不但可以提高计算精度,还可以有效避免沙漏变形.然而采用2×2×2高斯积分所引起计算量的增大是侵彻贯穿数值计算中必须面对的现实问题.大量冲击力学数值计算的实例表明,出现网格严重畸变,导致沙漏变形发展失控的单元往往只是部分单元,甚至是个别单元,因此我们可以采用局部2×2×2高斯积分,将单点积分计算效率高和全高斯积分可以消除沙漏变形的优点结合起来.在数值模拟侵彻贯穿问题中,对部分网格严重畸变单元采用2×2×2高斯积分,而其他的绝大部分单元仍采用单点积分,计算实践表明,这一方法是有效的.本文对八节点立方元分别采用不同积分方法对平头弹侵彻靶板进行了数值模拟计算,表明局部2×2×2高斯积分能够很好地控制沙漏变形,精度接近2×2×2高斯积分,且具有较高的计算效率.还采用单点积分和2×2×2高斯积分分别计算了板中一维应变波的传播,并与理论解进行了比较,结果表明2×2×2高斯积分的计算结果明显优于单点积分.1　计算方法简介1.1　单点积分的沙漏变形与沙漏控制目前,在大多数冲击力学有限元程序里,为提高计算效率,在计算节点力时,通常采用的是单点积分法,对于八节点六面体立方元(图1)而言,即在自然坐标系里,在进行单元体积分时,积分点取为ξ=η=ζ=0.单点积分的直接后果是立方元被认为是常应力常应变单元,因此其对应的位移场只能是线性位移场,显然,单点积分无法描述非线性位移模态.当单元变形出现如图2所示的几种沙漏变形模态时,由单点积分定义的单元应变为零,这就是所谓零能模式.在这种变形模式中,相应的节点内力 f int i= 0.于是沙漏变形可以在无节点阻力的情况下任意发展,在数值计算中它将导致网格相交,甚至网格翻转,计算中断.针对单点积分中出现的这一问题,Belyt schko 等[3,4]提出了一种通过引入抗沙漏节点力的方法控制沙漏变形.其基本思想是根据沙漏变形模态引入相应的沙漏变形率,q iα=18v iIγαI.(1)式中,角标i表示空间坐标轴方向,I表示节点,α表示沙漏变形模态;v iI为节点速度分量,γαI为沙漏变形分量,定义为γαI=ΓαI-1VB iI x iJΓαJ.(2)式中,B iJ=5N I(x i)/5x i,ΓαI如表1所示.由虚功原理出发,可以建立抗沙漏节点力,f H G iI=18Q iαγαI.(3)233中国科学技术大学学报第39卷式中,Q i α是与沙漏变形相对应的沙漏应力,Q i α=k λ+2μ3B iI B iIVq i α,(4)k 为人工沙漏控制系数,λ和μ是拉梅常数.1.2　2×2×2高斯积分的计算过程单点积分所对应的位移场是线性位移场,实际上八节点立方元所设定的位移函数应包含若干非线性项,或者说在自然坐标系里,形函数应具有如下的形式,N I (ξ,η,ζ)=18(1+ξξI )(1+ηηI )(1+ζζI )= 18{ΣI +ξΛ1I +ηΛ2I +ζΛ3I +ηζΓ1I + ζξΓ2I +ξηΓ3I +ξηζΓ4I }.(5)式中,ξI ,ηI ,ζI ,ΣI ,Λ1I ,Λ2I ,Λ3I ,Γ1I ,Γ2I ,Γ3I ,Γ4I 分别表示与节点相对应的自然坐标取值,如表1所示,其中,Γ1I ,Γ2I ,Γ3I ,Γ4I 刻画了位移场的非线性模态.表1　节点的自然坐标取值T ab.1　N atural coordinates of nodes节点ξηζΣIΛ1I Λ2I Λ3I Γ1IΓ2IΓ3IΓ4I1-1-111-1-11-1-11121-1111-11-11-1-13111111111114-1111-1111-1-1-15-1-1-11-1-1-1111-161-1-111-1-11-1-11711-1111-1-1-11-18-11-11-11-1-11-11采用有限元方法离散运动方程时,通常取空间域和时间域分别离散的所谓半离散化思想,空间域采用有限元方法离散,时间域则按有限差分方法离散.由半离散化思想知,速度场可表示为v ( x ,t )=N ～( x )・ u (t ).(6)式中,N ～( x )是由形函数构成的3×24矩阵,N ～( x )=N 100…N 8000N 10…0N 800N 1…N 8;u (t )表示由节点速度构成的24×1列向量,u (t )=[u 11　u 21　u 31　u 12　u 22　u 32　…　u 18　u 28　u 38]T,u ij 表示第j 个节点i 方向的速度分量. u (t )只是时间t 的函数.于是,v ・( x ,t )=N ～( x )・ u ・(t ).(7) 应变率矢量可以表示成ε・= ε11 ε22ε33γ12γ23 γ31=55x 100055x 2055x 355x 255x 1055x 355x 255x 355x 1v 1(x ,t )v 2(x ,t )v 3(x ,t )=L ・ v (x ,t ).(8)式中, εii 为法向应变率, εii =5v i5x i; γij 为剪应变率, γij =5v i 5x j +5v j5x i; L 是由微分算子构成的矩阵,v i (x ,t )表示x i 方向上的速度分量.将式(6)代入式(8),有ε・= L ・N ～( x )・ u (t )= B ( x )・ u (t ).(9)式中, B 称为梯度矩阵,由形函数的空间导数构成,B ( x )= L ・N ～( x ).由式(5)可知,形函数N i (ξ,η,ζ)是由自然坐标给出的,根据偏微分法则,5N i 5ξ=5N i 5x 15x 15ξ+5N i 5x 25x 25ξ+5N i 5x 35x 35ξ.(10)对于其余两坐标(η,ζ),可以类似得到,写成矩阵形式有5N i5ξ5N i5η5N i 5ζ=5x 15ξ5x 25ξ5x 35ξ5x 15η5x 25η5x 35η5x 15ζ5x 25ζ5x 35ζ5N i5x 15N i 5x 25N i 5x 3=J 5N i 5x 15N i5x 25N i 5x 3(11)式中,J 称为雅可比矩阵,可记为5(x 1,x 2,x 3)5(ξ,η,ζ).在等参元里,空间坐标可以表示为x j =∑N ixji=N 1x j 1+N 2x j 2+…+N 8x j 8;i =1,2,…,8;j =1,2,3.(12)式中,x ji 表示第i 个节点j 方向的坐标值.利用式(12),J 可以显式地表现为自然坐标的函数,J ≡5(x 1,x 2,x 3)5(ξ,η,ζ)=333第3期冲击力学有限元计算中的数值积分∑5N i5ξx 1i ∑5N i5ξx 2i ∑5N i5ξx3i ∑5N i5ηx 1i ∑5N i5ηx 2i ∑5N i5ηx 3i ∑5N i5ζx 1i ∑5N i5ζx 2i ∑5N i5ζx 3i=5N 15ξ5N 25ξ…5N 85ξ5N 15η5N 25η…5N 85η5N 15ζ5N 25ζ…5N 85ζx 11x 21x 31x 12x 22x 32………x 18x 28x 38.(13)于是N i 对于(x 1,x 2,x 3)的偏导数可用自然坐标表示为5N i 5x 15N i 5x 25N i5x 3=J -15N i 5ξ5N i5η5N i5ζ.(14)式中,J -1是J 的逆矩阵.体积微元可以表示为d V =d ξ・(d η×d ζ).(15)式中,d ξ=5x 15ξd ξ i +5x 25ξd ξ j +5x 35ξd ξk ,d η=5x 1ηd η i +5x 2ηd η j +5x 3ηd η k ,d ζ=5x 15ζd ζ i +5x 25ζd ζ j +5x 35ζd ζ k ; i , j , k 是直角坐标的单位向量,于是d V =5x 15ξ5x 25ξ5x 35ξ5x 1η5x 2η5x 3η5x 15ζ5x 25ζ5x 35ζd ξd ηd ζ=|J |d ξd ηd ζ.(16) 节点内力可以表示为f inti (t )=∫ΩeB T i ( x ) σ( x ,t )d Ω= ∫VB T i (ξ,η,ζ) σ(ξ,η,ζ,t )d V = |J |∫1-1∫1-1∫1-1B Ti(ξ,η,ζ) σ(ξ,η,ζ,t )d ξd ηd ζ.(17)采用2×2×2高斯积分[5]后,式(17)可写为 f int i (t )=|J |∑2j =1∑2k =1∑2m =1HjH k H m B T i (ξj ,ηk ,ζm ) σ(ξj ,ηk ,ζm ,t ).(18)式中,B Ti σ=5N i5x 1005N i5x 205N i 5x 305N i5x 205N i 5x 15N i 5x 305N i 5x 305N i 5x 25N i 5x 1σ11σ22σ33σ12σ23σ31=5N i 5x 1σ11+5N i 5x 2σ12+5N i 5x 3σ315N i 5x 2σ22+5N i 5x 1σ12+5N i 5x 3σ235N i 5x 3σ33+5N i 5x 2σ23+5N i 5x 1σ31.(19)对于2×2×2高斯积分,H j =H k =H m =1,ξj =ηk =ζm =±0157735(j ,k ,m =1,2),将式(14)和式(18)代入式(17)可以计算得到节点内力.由于在每个坐标方向的积分点数取为n =2,数学上可以证明,2×2×2高斯积分的结果可以达到2×n -1=3阶的精度.也就是说,如果 B T i (ξ,η,ζ) σ(ξ,η,ζ)=∑ξi ηj ζk,且i ,j ,k ≤3,则式(18)将能给出积分的精确解.2　数值模拟计算及结果分析2.1　平头杆弹对靶板的三维侵彻计算为了考查沙漏变形对有限元计算过程的影响,我们对平头杆弹侵彻靶板进行了数值计算.弹体和靶板材料取为线性硬化体,其静水压力与体积变形之间用如下的Gruneisen 状态方程描述[6]:P =(K 1μ+K 2μ2+K 3μ3)(1-0.5γμ)+γρE.(20)式中,μ=V 0/V -1为压缩比,γ为Gruneisen 参数,E 为比内能.应力偏量服从增量型弹塑性应力应变关系,即d S ij =2G (de ij -d e p ij).(21)材料服从Mises 屈服准则,屈服强度取为线性硬化函数,Y =Y 0+(Y u -Y 0) εP / εu .(22)式中,Y 0为初始屈服应力,Y u 为极限屈服应力, εP433中国科学技术大学学报第39卷表2　材料参数T ab.2　P arameters of m aterialρ/(kg ・m -3)G /GPaY 0/GPaY u /GPaP σ/GPa εu状态方程参数/GPaK 1K 2K 3γ782377.5 1.43 2.50 1.700.011642945002.图3　不同积分方法下的侵彻图像Fig.3　Penetration conf igurations of different integration methods为等效塑性应变, εu 为极限应变.计算中用到的有关物理参数[6]和计算参数如表2所示.图3给出了平头杆弹以初速度V 0=1500m/s 侵彻靶板t =2710μs 时的应变及网格变形图像.可以看出,加入抗沙漏节点力的单点积分和2×2×2高斯积分都能够很好地控制沙漏变形,避免网格畸变,而取消抗沙漏节点力的单点积分则出现了网格的畸变和相交,严重影响了有限元计算的正常计算过程,这也进一步说明了上述计算方法的正确性.尽管2×2×2高斯积分的计算结果可以更为准确地给出杆弹靶板的物理图像,然而其计算量的增加依然是数值模拟高速侵彻和贯穿问题时不可忽视的因素.事实上,计算中我们发现,出现网格畸变严重导致沙漏变形模态发展的只是部分单元,这部分单元几乎全部集中在侵彻通道附近.为此我们开展了一组只对侵彻通道附近靶板单元采用2×2×2高斯积分,其余单元都采取不加抗沙漏节点力的单点积分,结果如图3(d )所示.由于充分利用了单点积分计算效率高、2×2×2高斯积分可以比较准确捕捉网格变形的物理细节等两者的长处,因此这种局部2×2×2高斯积分方法是非常有效的.图4是平头杆弹以初速度V 0=1500m/s 侵彻靶板时的平均侵彻速度与时间的关系曲线.可以看到,加入抗沙漏节点力的单点积分、2×2×2高斯积分和局部2×2×2高斯积分的模拟结果基本一致,而没有抗沙漏节点力的单点积分得到的平均侵彻速度在相同时刻远大于前两者,并且存在局部震荡,这是由于后者靶板单元产生过大的网格畸变而导致破图4　不同积分方法下的平均侵彻速度随时间的变化Fig.4　History of average velocity for different integration methods坏严重的缘故.表3给出了单点积分、2×2×2高斯积分和局部2×2×2高斯积分在不同的条件下所用CPU 机时的比较,可见与2×2×2高斯积分相比,局部2×2×2高斯积分的效率较高,且所用的计算时间只是稍大于单点积分.我们认为,以稍多的计算时间来换取更高的计算精度是值得的.表3　不同积分方法所用CPU 机时比较(单位:秒)T ab.3　CPU time used by different integration methods(unit :second)方法三维侵彻计算的靶板网格数12003000720010000单点积分5321896567696252×2×2高斯积分73531221073226354局部2×2×2高斯积分5952056601510463533第3期冲击力学有限元计算中的数值积分图5　脉冲载荷下单点积分和2×2×2高斯积分的应力波形Fig.5　Stress w aveform of single point qu adrature and2×2×2G aussian qu adrature in pulse load2.2　一维应变波的对比计算为确认上述积分方法的有效性和合理性,我们采用在此基础上开发的三维Lagrange有限元程序计算了脉冲载荷下半无限板中的弹塑性波的传播问题.设厚L=014m的半无限厚钢质平板上作用有脉冲载荷P,P=P0,t≤τ;0,t>τ.(23)式中,P0=5×109Pa,τ=10μs.当P0大于材料的屈服强度时,平板内部将产生弹塑性间断波的传播.平板材料的本构模型和相关参数同节211.图5给出了两种积分方法所得的计算结果和解析解的比较.当加载应力大于材料的初始屈服强度时将形成双波结构,较慢的塑性激波尾随在较快的弹性前驱波之后.随着时间的推移,弹性前驱波与塑性主波阵面逐渐拉开,而卸载时产生的弹性卸载波又对塑性波产生追赶卸载.由图5可见,单点积分的精度较低,尤其是在间断面上,其波形失真较大,间断面被过度地光滑化了.和单点积分相比,2×2×2高斯积分的计算精度有较大的改善,得到的波形更接近于理论波形.因此,通过数值计算方法模拟应力波演化过程时,最好放弃单点积分而采用全高斯积分方式.3　结论本文在Lagrange有限元数值计算方法的基础上,介绍了两种用于计算节点力的不同积分方法,并详细阐述了2×2×2高斯积分法的计算过程.单点积分计算效率高,但需要进行沙漏控制;2×2×2高斯积分可以避免沙漏变形,且计算精度较高,但计算量大;而采用局部2×2×2高斯积分则同时具有两者的优点.采用自建的三维有限元程序,对平头弹侵彻靶板进行了数值模拟计算,表明局部2×2×2高斯积分能够很好地控制沙漏变形,精度接近2×2×2高斯积分,且具有较高的计算效率.板中一维应变波传播的数值模拟对比计算也表明2×2×2高斯积分的计算结果明显优于单点积分.进一步的研究可以尝试采用其他的单元类型,例如四节点四面体单元,来进行网格的剖分,这样可以从根本上避免沙漏变形.当然,这会引起滑移面计算、积分方法等一系列计算方法的改变.参考文献(R eferences)[1]Znkas J A.High Velocity Impact Dynamics[M].NewY ork:Wiley,1990.[2]Camacho G T,Ortiz M.Adaptive Lagrangian modelingof ballistic penetration of metallic targets[J].Comput Meth Appl Mech Eng,1997,142(324):2692301.[3]Flanagan D P,Belytschko T.A uniform strainhexahedron and quadrilateral with orthogonal hourglass control[J].International Journal for Numerical Methods in Engineering,1981,17:6792706.[4]Belytschko T.Correction of article by D P Flanaganand T Belytschko[J].International Journal for Numerical Methods in Engineering,1983,19:4672468.[5]王勖成.有限单元法[M].北京:清华大学出版社,2003.[6]Robbins J R,Ding J L,Gupta Y M.Load spreadingand penetration resistance of layered structures:A numerical study[J].Int J Impact Eng,2004,30(6): 5932615.633中国科学技术大学学报第39卷。

1．针对下图所示的3个三角形元，写出用完整多项式描述的位移模式表达式。

2．如下图所示，求下列情况的带宽:a)4结点四边形元；b)2结点线性杆元。

3．对上题图诸结点制定一种结点编号的方法，使所得带宽更小。

图左下角的四边形在两种不同编号方式下，单元的带宽分别是多大？4．下图所示，若单元是2结点线性杆单元，勾画出组装总刚后总刚空间轮廓线。

系统的带宽是多大？按一右一左重新编号（即6变成3等）后，重复以上运算。

5． 设杆件1－2受轴向力作用，截面积为A ，长度为L ，弹性模量为E ，试写出杆端力F 1，F 2与杆端位移21,u u 之间的关系式，并求出杆件的单元刚度矩阵)(][e k6．设阶梯形杆件由两个等截面杆件○1与○2所组成，试写出三个结点1、2、3的结点轴向力F 1，F 2，F 3与结点轴向位移321,,u u u 之间的整体刚度矩阵[K]。

7． 在上题的阶梯形杆件中，设结点3为固定端，结点1作用轴向载荷F 1=P ，求各结点的轴向位移和各杆的轴力。

8． 下图所示为平面桁架中的任一单元，y x ,为局部坐标系，x ，y 为总体坐标系，x 轴与x 轴的夹角为 。

（1） 求在局部坐标系中的单元刚度矩阵 )(][e k （2） 求单元的坐标转换矩阵 [T]；（3） 求在总体坐标系中的单元刚度矩阵 )(][e k9．如图所示一个直角三角形桁架，已知27/103cm N E ⨯=，两个直角边长度cm l 100=，各杆截面面积210cm A =，求整体刚度矩阵[K]。

10． 设上题中的桁架的支承情况和载荷情况如下图所示，按有限元素法求出各结点的位移与各杆的内力。

11． 进行结点编号时，如果把所有固定端处的结点编在最后，那么在引入边界条件时是否会更简便些？12． 针对下图所示的3结点三角形单元，同一网格的两种不同的编号方式，单元的带宽分别是多大？13． 下图所示一个矩形单元，边长分别为2a 与2b ，坐标原点取在单元中心。