拉格朗日插值公式

2024-07-15

拉格朗日插值公式（精选3篇）

拉格朗日插值公式篇1

现代电力系统中,架空输电线路扮演着重要角色,它不仅承担着传送输配电的功能,还要抵抗自然或人为带来的干预与破坏,且需要对铁塔的运行状态给出合理的综合评价,防患于未然,提高输电线路运行的安全性和可靠性。长期运行资料表明,输电线路安全事故大多由输电铁塔本体受外力破坏造成[1]。输电铁塔本体安全由铁塔应力和塔材实际强度2个方面表征[2]。随着我国电力事业的迅速发展,输电铁塔结构日趋复杂,对于存在柔性杆件的铁塔,若采用传统的有限元线性分析,则计算结果由于计算过程中对其进行线性化近似处理而产生较大的偏差。目前采用有限元非线性分析[3,4,5,6],该方法将铁塔结构中的刚性杆件也用非线性分析,虽然其计算结果的精确度较高,但大大增加了对刚性杆件分析的计算量,降低了计算速度,不能应用于铁塔结构安全的实时监测。因此本文将铁塔结构中的刚性单元跟柔性单元分别进行线性分析与非线性分析,以提高其计算速度。

目前对铁塔结构的应力分析都采用有限元法,有限元是将所研究的物体分解成若干个单元,每个单元先假定一个近似解,然后求出该域的满足条件,进一步得到最终的解[7,8]。有限元与传统的分析方法相比具有较高的精度,且能分析复杂的铁塔结构,然而有限元只能求解节点的应力,并通常把节点应力视为铁塔构件的最大应力,实际情况中构件的最大应力并非就是节点应力,且构件的最大应力由该构件全部节点变形计算得到,而节点应力只由该节点变形计算得到,因此若将节点应力视为最大应力则会影响铁塔结构安全评价的准确性。本文对铁塔结构中的刚性单元跟柔性单元分别采用线性和非线性方法对其进行分析,且对节点应力进一步分析处理,求出2节点间构件的最大应力,得到铁塔应力的精确分布及最大值,提高对铁塔结构安全评价的准确性。

1 铁塔结构模型与材料模型

1.1 铁塔结构模型

首先建立铁塔结构的桁梁混合模型,将既承受轴向力又承受剪力和弯矩的主材或者横隔材视为梁单元,只承受轴向力的斜材被视为杆单元,而不承受作用力的辅材则被简化掉不作为模型的单元[9,10,11]。以铁塔的横担方向作为整体坐标系的x轴,线路方向作为y轴,竖直方向作为z轴,并满足右手定则;以杆单元所在直线作为单元局部坐标系的x轴,杆件与局部坐标系下的x轴方向重合,其正方向与整体坐标系x轴正方向一致。

1.2 铁塔材料模型

对于不存在柔性杆件的铁塔,塔材是线弹性材料。对于存在柔性杆件的复杂铁塔,将塔材分为2组:(1)承受拉压的刚性单元:塔材是线弹性材料。(2)承受拉力的柔性单元:只承受拉力,不能承受压力杆件[12,13]。

2 铁塔结构应力分析

2.1 刚性单元线性分析

根据所用钢材的横截面面积、材料的弹性模量和剪切模量等相关参数生成单元刚度矩阵[14],根据杆件之间的空间角度关系、杆件之间连接关系,转换叠加出铁塔整体刚度矩阵,根据铁塔钢材的自重、铁塔所受的风荷载和铁塔所受导线的拉力分别等效为相应节点所受的载荷,并生成节点载荷阵列,每个节点具有6个自由度,即杆件承受一维轴力、两维剪力、两维弯矩、一维扭矩,即对应着节点的6个自由度[14];以节点位移阵列作为未知量,与整体刚度矩阵,节点载荷阵列组成矩阵方程。由于整体刚度矩阵为奇异矩阵,方程组无解,若要求解该方程,必须引入约束条件,限制铁塔结构的刚性位移,保证整体刚度方程有惟一解。采用对角元素置1法,将δi=δ0引入整体刚度矩阵,针对输电铁塔的4个塔腿中,与基础连接的部分是固定端约束,因此δ0=0;将刚度矩阵K的第i行的主对角线元素Kii置1,其余元素清0,且将第i行的载荷项Ri用0代替。即代入24个位移边界条件,消除整体刚度矩阵的奇异性,从而采用高斯消元法进行矩阵方程求解,求解出节点位移矩阵,再根据弹性力学中应变与位移的关系[15],计算出各节点的应力、应变。

2.2 柔性单元非线性分析

由于铁塔结构中柔性杆件的应力与应变呈非线性关系。因此,对于求解此类非线性问题,不能采用传统的直接求解方法,必须把非线性问题分成若干个加载步,分阶段对其逐步求解,即只要把荷载分的足够细,迭代次数足够多,就可以用分段线性分析代替大位移小应变的非线性[16]。运用修正的结构几何位置变形原理对铁塔结构中的柔性单元进行非线性分析,即以t时刻的状态作为基准,推出t+Δt时刻的状态。修正的结构几何位置变形原理:

式(1)中:X0,Y0为柔性单元的坐标值;Ui,Vi为单元变形后在节点i处的位移;Uj,Vj为单元变形后在节点j处的位移。

梁单元的节点位移可以表示为:

式(2)中:l0,θ0为柔性单元的坐标值;Ui,Vi,θi为单元变形后在节点i处的位移;Uj,Vj,θj为单元变形后在节点j处的位移。

则节点位移阵列可以表示为:

变形后的单元节点力可以表示为:

式(4)中:[K'](e),[δ'](e),[F'](e)分别为局部坐标系下变形后的单元刚度矩阵、节点位移和单元节点力。

通过坐标转换为整体坐标下单元节点力,式(4)可变为:

如果将结构以线性分析计算得到的弹性位移作为第一次近似值,然后通过式(3)、式(4)算出各单元作用在节点上的力为:

则在各节点上产生的不平衡力为:

将不平衡力作用到结构的各节点上,得出节点的第二次近似值,重复上述过程多次迭代直至[ΔR]≈0为止。假设结构在载荷作用下已用线性理论方法求出位移的近似值,其迭代步骤为:

(1)建立各单元的局部坐标,并计算出各单元在局部坐标下的单元刚度矩阵[K'](e)和位移阵列[δ'](e)。

(2)计算节点应力,将局部坐标下的单元刚度矩阵[K'](e)和节点应力[F'](e)经坐标变换转换到整体坐标下的

(3)叠加生成整体结构的刚度矩阵

(4)计算出各单元作用于节点上的力[Rr],并计算不平衡力[ΔR],即:

(5)重复上述过程多次迭代直至[ΔR]≈0为止。

2.3 拉格朗日插值

根据上述计算方法所得到的节点应力、应变对各个矩阵中的各项值进行拉格朗日插值,通过插值函数的计算得到较为精确的铁塔各杆件的应力计算公式。应力的拉格朗日插值表达式为:

式(9—14)中:Fi为节点的应力矢量;x,y,z分别为节点应力的方向;li为拉格朗日基本多项式(拉格朗日基函数);L为拉格朗日插值多项式。

对其中的L(x),L(y),L(z),L(xy),L(yz),L(xz)的自变量进行一阶微分,求出其导数等于0的点,即令L'(x)=0,L'(y)=0,L'(z)=0,L'(xy)=0,L'(yz)=0,L'(xz)=0,解分别记为x',y',z',xy',yz',xz',分别求出L(x'),L(y'),L(z'),L(xy'),L(yz'),L(xz')的值,此时可以求得铁塔x,y,z轴各个方向上应力的极点以及最大值。

3 设计与实现

以桁梁混合模型对铁塔结构建模,将既承受轴向力又承受剪力和弯矩的主材或者横隔材视为梁单元,只承受轴向力的斜材被视为杆单元,而不承受作用力的辅材则被简化掉不作为模型的单元。对于存在柔性杆件的复杂铁塔,将铁塔结构按材料分为2组:刚性单元和柔性单元,对于刚性单元采用有限元线性分析,即上述求解矩阵方程的方法求出节点的应力、应变;对于柔性单元采用有限元非线性分析,即利用上述的有限元增量法和牛顿迭代法求解节点的应力与应变。根据计算得到的节点应力、应变,分别对其进行拉格朗日插值,通过插值函数的计算得到较为精确的铁塔整体各部位的应力计算公式,并对其进行一阶微分求出最大值,即铁塔构件中的最大应力。在eclipse的开发环境下,用java语言编写铁塔构件应力精确计算程序,从而实现上述的计算方法,提高对输电铁塔结构安全评价准确性,计算流程图如图1所示。

4 实例验证

110 k V直线塔倒塔事故照片如图2所示。

以图2中的110 k V直线塔为例,利用上述的计算方法,对直线塔结构中的刚性单元进行线性分析,对柔性单元进行非线性分析,求出各节点的应力,并对其进行拉格朗日插值,求出铁塔构件中的最大应力,把计算出的构件最大应力与把节点应力视为铁塔构件的最大应力以及钢材的屈服强度进行比较,如表1所示。表1左侧最大应力为本方案计算结果,右侧最大应力则为将节点应力视为铁塔构件的最大应力的计算结果。

由表1可知,传统有限元分析中,塔身部分编号为5-12,6-10,7-11,11-12,9-12,14-15,15-16,13-16的8根杆件的应力未超过其屈服强度,然而在本方案计算结果中却超过了屈服强度,因而导致杆件发生弯曲变形、折断,造成钢材的变形折断,与图2实际倒塔事故照片中铁塔折断的位置吻合,计算结果更加精确,因此可以利用该计算方法对输电铁塔结构进行更准确的安全评价。

5 结束语

本文从理论和实际工作2个方面对输电铁塔构件的最大应力分析做了初步的探讨,基于铁塔结构的有限元原理分析方法,自主开发了一种计算输电铁塔构件最大应力的方法,并开发了输电铁塔应力精确计算的软件。通过上述的实例验证,证明了该方法的计算结果更为精确,可以提高对铁塔结构安全评价的准确性,并可以对铁塔最薄弱的环节进行预警。

拉格朗日插值公式篇2

地形形态特征分析主要是在地面高程数据中提取相关的地形特征信息,这些地形特征往往受地面高程数据数据结构的影响不能在相应的数字高程模型(DEM)中直接表达,其中以格网DEM的影响最为明显[1]。加之近年来自然灾害频繁发生,地形形态特征分析受到越来越多的关注,流域分析作为地形形态特征分析所涉及的一个研究课题也受到了空前的关注,流域分析的基础和重点即地形形态特征分析,从流域分析角度看,地形形态特征分析的核心即地形形态特征线提取,从算法原理上讲,地形形态特征线提取包括两类基本方法:解析法和模拟法。模拟法主要根据地表物质水流运动特性,利用水流模拟提取地形特征线(即分水线与汇水线),流域分析中地形特征线的提取正是利用了流域汇水面积的特性实现了分水线与汇水线的提取,在DEM上进行汇水面积计算的算法称为水流路径算法,水流路径算法的实现与DEM结构有关。前文已提到,格网DEM对地形形态特征的影响最为明显,因而本文重点对格网结构DEM进行水流路径分析。迄今为止,基于格网DEM已出现众多水流路径算法,由于众学者对水流向下游的流量分配是否单一的看法不一,形成单流向和多流向两类水流路径算法,由于单流向算法计算相对简单、效率较高且对凹地、平坦区域有较强的处理能力而使其应用较为广泛,具有代表性的单流向算法有D8算法[2],随机方向算法[3]和流向驱动算法等。在实际地形形态特征线提取中,如土丘中的细长冲沟等处,由于DEM分辨率和水流路径算法模拟精度等的限制,致使推导的水流方向可能会产生较大误差,即当前点应流入下游唯一的高程最低格网点,此时若采用随机方向算法不仅会直接影响水流方向和局部区域汇水面积的精确计算,甚至导致最终提取的分水线和汇水线在局部区域内也会产生较大误差。为解决此类问题,提出了一种水流路径算法,该算法通过构建经过当前格网点的局部曲面计算该点在局部曲面中的等高线曲率,并利用所求得的等高线曲率选取合适的单流向算法来推断当前点的水流流向,进而计算出整个区域的水流累计量。

1 局部曲面拟合与等高线曲率的计算

局部曲面拟合是DEM中常用到的方法,它是一个连续且光滑的数学曲面,局部曲面拟合是根据空间抽样数据拟合一个数学曲面,用该数学曲面反映空间分布的变化情况。目前已有许多成熟的局部曲面拟合函数,如二元样条函数、多项式函数、Geomap曲面、多层曲面叠加、最小二乘配置方法等[4]。由于曲面由曲线按一定规律和法则交叉排列构成,且曲面拟合的基础即曲面点和线的规则组合,所以本文首先根据局部格网点构建局部曲线,进而利用局部曲线拟合局部曲面。

1.1 局部曲面拟合

基于格网DEM进行流域分析的基础即区域格网点水流方向的确定,水流方向是水流离开格网流向其下游格网点时的方向。在流域分析中,通常在3×3的局部窗口中按照一定的水流路径算法确定水流方向,此处通过在局部窗口中构建局部曲面并对其进行分析实现水流方向的确定,局部曲面函数通过拉格朗日公式构建的6条曲线按一定权值组合而成,在3×3局部窗口中,各行各列都可组成一条拉格朗日曲线,每行基于X值变化(即i)和高程Z可构成一条拉格朗日曲线,每列基于Y值变化(即j)和高程Z也可构成一条拉格朗日曲线,3条基于X值的曲线根据一定权值可构成描述局部曲面X值的局部曲面子函数曲线,3条基于Y值的曲线也可根据权值构成描述局部曲面Y值的子函数曲线,最后两个拉格朗日子函数曲线按一定权值构成局部窗口内的局部曲面,如图1。

第一、二、三行3个格网点X值与Z值构成的拉格朗日曲线为:

第一、二、三列3个格网点Y值与Z值构成的拉格朗日曲线为:

3条基于X值和Z值的拉格朗日曲线在给定权值因子Pxz1、Pxz2和Pxz3的情况下可构成描述局部曲面X值和Z值的局部曲面子函数曲线LXZ(x)。同理3条基于Y值和Z值的拉格朗日曲线根据权值因子也可构成描述局部曲面Y值和Z值的子函数曲线LYZ(y),如下:

局部曲面子函数LXZ(x)和LYZ(x)求出后,引入两个影响地形拟合精度的权值因子PX和PY,并根据其构建出最终的局部曲面函数Z:

其中PX=max(Pxz1,Pxz2,Pxz3),PY=max(Pyz1,Pyz2,Pyz3),即PX和PY要选择最符合区域地形的行曲线和列曲线的XZ和YZ权值因子来对其进行代替。

1.2 计算等高线曲率

局部曲面中的等高线曲率就是通过该点的等值面(水平面)与地表交线的曲率(即通过该点的等高线的曲率),一定程度上等高线曲率表达了地表物质(如水流等)运动的聚合和发散程度[5],笔者认为等高线曲率既然能较好的描述水流聚合与发散状况,那么可基于此来模拟推断当前点的水流方向,通过局部曲面某点的等高线曲率计算公式如下:

2 水流路径算法研究

在实际地形形态特征分析中,尤其是对流域分析中的地形特征线进行提取时,在局部地形的某些区域如土丘中的细长冲沟、丘陵中的侵蚀细沟以及微地貌中的侵蚀犁沟和小河沟的细小河谷等地,常常由于DEM分辨率的限制或水流流量阀值的不当造成水流方向在这些特殊地形处的误判,并直接影响到局部区域汇水面积的精确计算,进而导致提取的分水线和汇水线也会产生一定误差,如当前点在细长冲沟处的水流方向为该点与冲沟中与该点临近的下游唯一的高程最低点的方向,很明显,此时若选用随机方向算法就很有可能造成水流方向在这些特殊地形特征点处的误判。为减少此类问题的出现,提出了一种水流路径算法,该算法以当前格网点为中心,建立一个3×3的局部窗口,然后提取窗口内对应格网点的几何信息,根据这些已知数据构建出局部曲面函数并模拟推断该点的水流方向,具体步骤如下:

2.1 以当前格网点(i,j)为中心,利用3×3局部窗口中提取的9个格网点的(X,Y,Z)坐标根据式(1)-(4)构建出基于该局部窗口的局部曲面函数Z=f(x,y),再根据式(5)求出当前点在局部曲面中的等高线曲率CC。

2.2 设置一个等高线曲率阀值C(设C>0):

(1)若CC∈(-C,C),证明当前点的水流聚合程度较高,发散程度较低,则当前点的水流方向应该唯一的指向临近该点的下游最低高程的格网点,即此时采用D8算法确定当前格网点的水流方向及流量分配比例;

(2)若CC∈(-∞,-C]∪[C,+∞),证明当前点水流聚合程度较低,而水流发散程度相对较高,促使该点水流不一定流入下游最低高程格网点,则此时采用随机八方向算法(Rho8)确定当前格网点的水流流向及其向下游的流量分配比例。

当前格网点的水流方向及流量分配比例确定后,再按一定搜索法则搜寻局部区域格网DEM中的下一格网点,依次循环进行,便可计算出整个DEM研究区域的水流方向矩阵,进而可计算出研究区域内的汇水面积,并为区域分水线和汇水线的成功提取奠定基础,图2和图3为C=18.6时提取的水流方向和水流累积量。

3 结论

以当前格网点为中心构建局部曲面,通过该点在局部曲面中的等高线曲率可选取合适的单流向算法来推断当前点的水流方向,这为区域分水线和汇水线的提取奠定了基础。该方法在流域分析中已取得较好效果,提取的分水线和汇水线也符合精度要求,但由于要根据等高线曲率的动态变化来选取不同的单流向算法,致使对格网DEM进行数据处理时速度较慢、效率较低,需进一步研究,此外等高线曲率阀值C也需根据实际地形趋势适当选取。

摘要：提出了一种水流路径算法,该算法通过以各格网点为中心建立3×3窗口内基于6条拉格朗日曲线构成的局部曲面,再利用当前点在局部曲面中的等高线曲率选取合适的单流向算法以确定当前点的水流方向及其向下游点的流量分配比例。最后通过实例证明了该方法在水流累积量计算中的良好效果。

关键词：DEM,单流向算法,拉格朗日公式,等高线曲率

参考文献

[1]周启鸣,刘学军.数字地形分析[M].北京:科学出版社,2006.

[2]O’Callaghan J F,Mark D M.The Extraction of Drainage Networks from Digital Elevation Data[J].Computer Vision,Graphics,and Image Processing,1984,28:323-344.

[3]Fairfield J,Leymarie P.Drainage Networks from Grid Elevation Models[J].Water Resources Research,1991,27(5):709-717.

[4]汤国安,刘学军,闾国年.数字高程模型及地学分析的原理与方法[M].北京:科学出版社,2005.