1.2 工程问题的数学物理方程
下面通过求解微分方程的变分形式,给出微分方程的解。由于物理问题的边界条件是由变分形式自然给出的,因此求解问题的变分形式常常包含了某些边界条件。
为了便于阐述基本概念,所论述的变分函数将不考虑边界条件,典型微分方程的边界条件由相应的物理条件给出。
1.2.1 工程问题的数学物理方程概述
工程问题的控制方程通常由其基本方程和平衡方程给出。下面具体介绍诸多不同物理问题的一维控制方程。控制方程在基本形式上是大体一致的,并采用与相关工程领域中的用法一致的数学符号。
从变分原理出发,可以给出相应的物理问题的有限元模型。对于质量、力、长度、时间、温度及能量等物理量,分别使用M、F、L、t、T和E来表示。
1.一维弹性问题
在法向应力和轴向外力的作用下,直杆上的力的平衡问题可以表示为一维微分方程问题。假设x处的应力分布为σ(x),截面积为A(x),轴向外力为f(x),则直杆在该处的应力为σ(x)A(x),考虑到外力对直杆的作用,可以得到如下的基本方程:
由胡克定律(Hooke)可知,应力和截面x处的伸长率(即应变)ε之间的关系由构成直杆材料的弹性系数ε(x)给出。利用弹性模量(杨氏模量)定律,应力与轴向位移μ(x)之间满足如下关系:
或者
联合式(1-2)和式(1-3)可以得到如下关于位移μ(x)的二阶微分方程:
上述方程有两类不同形式的边界条件:一类为自然边界条件;另一类为几何边界条件,或者称为本质边界条件。其中,关于μ(x)的边界条件为本质边界条件;关于σ(x)的边界条件为自然边界条件。
相关物理量的单位:σ(x)为F/L2,A(x)为L2,f(x)为F/L3,ε为L/L,E(x)为F/L2,μ(x)为L。
将力的守恒方程与定义悬索的斜率的基本几何方程相结合,可以将发生微小弯曲后的悬索的平衡问题视为弹性问题。利用小扰动原理,假设悬索的张力T为常量,则作用在悬索上的合力满足如下方程:
式中,θ(x)为悬索的斜率,v(x)为悬索在垂直方向上的扰动量,k(x)为悬索的弹性模量,f(x)为作用在悬索上的、沿着垂直方向的载荷。利用小扰动原理,得知θ为一个极小量。从而可以得出悬索的基本几何方程,即θ可以被近似地表示为
由式(1-5)可以得出,张力沿着垂直方向的分量为
联合式(1-5)和式(1-6)可以得到如下的控制方程:
这里关于v(x)的边界条件为本质边界条件;而由式(1-7)给出的关于Fy的边界条件为自然边界条件。值得注意的是,上述边界条件应与给定的斜率θ等价。
相关物理量的单位:T为F,θ(x)为L/L,v(x)为L,f(x)为F/L,k(x)为F/L2。
2.热传导问题
从能量守恒方程和基本方程可以推导出描述一维定常热传导问题的基本方程。能量守恒原理要求热通量q的变化与外部的热源Q相等,即
式中,A(x)为传热面积,Q的相反数表示系统传出去的热量,其基本方程(也称为傅里叶定律)为
式中,T表示温度,k代表热传导系数。联合式(1-9)和式(1-10)得到如下的二阶控制微分方程:
式中,关于T的边界条件为本质边界条件;而关于q(x)的边界条件为自然边界条件。
相关物理量的单位:q(x)为E/tL2,A(x)为L2,Q(x)为E/tL3,T(x)为T,k(x)为E/tLT。
3.位势流
作为流体力学领域中的一个特殊方面,位势流广泛应用于地下水流动问题中,在上述应用中,可以假设流体为定常不可压缩流动,从而可以完全由连续方程或质量守恒方程描述上述问题。假设面积为一个定常数,则一元势函数可以假设为
或者
式中,μ为流动速度,h为压力头,z为高度头,γ为地下水问题中水的比重,p为压力,K为渗透率或水力学传导系数。其基本方程由达西定律(Darcy)给出:
上式表明达西定律与由式(1-12)给出的位势流的定义是相关的。注意到一维定常不可压缩流动满足du/dx=0,并联合式(1-12)~式(1-14)有
式(1-15)的解为线性函数,因而一维流动问题的速度为一个常数。然而,与之相应的二维问题的解却比较复杂。在上述问题中,关于φ的边界条件为本质边界条件;关于速度的边界条件为自然边界条件。
相关物理量的单位:φ(x)为L2/t,h(x)为L,K(x)为L/t,u(x)为L/t。
4.质量输运方程
对于大多数的基本位势流问题,如果其流动是定常的且控制方程与式(1-15)相似,则会发生扩散现象。在此类问题中,质量方程的平衡将以稀释混合项的形式写出。上述理论被广泛应用于物理问题中。特别地,当把地下水流动问题视为位势流问题时,可以假设某种物质与地下水构成混合物来共同流动,并同时在混合物中进行扩散,从而可以将位势流理论与质量扩散理论相结合,得到主要物理问题的一个较为全面的描述。假设截面积为一个定常数,则稀释混合项的质量可以表示为
式中,μ(x)表示混合物的流动速度;C(x)、j(x)分别表示稀释项的浓缩系数和流通量;Kr表示稀释项与其周围物质之间的反应速度,如化学反应;m为外部的质量源函数。其基本方程称为Fick定律,可以表示为
式中,D(x)为扩散系数,联合式(1-16)和式(1-17)可以得到如下的控制方程:
这里假设速度μ(x)是已知的。其中,关于C(x)的边界条件为本质边界条件;关于流通量j(x)的边界条件为自然边界条件。
相关物理量的单位:C(x)为M/L3,D(x)为L2/t,u(x)为L/t,Kr为t-1,j(x)为M/L2t。
5.电流
静电控制方程与热传导方程是相似的。下面的电荷平衡方程给出了电量分布D(x)与电荷密度ρ(x)之间的关系:
式中,A(x)为垂直于X轴的横截面面积,电场E(x)与电势φ(x)之间的关系满足:
相应的基本方程为
式中,ε(x)表示材料的电容。联合式(1-19)和式(1-21)得到如下的控制方程:
其中,关于φ的边界条件为本质边界条件;关于D的边界条件为自然边界条件。
相关物理量的单位:D(x)为Q/L2,A(x)为L2,ρ(x)为Q/L3,E(x)为V/L,φ(x)为V,ε(x)为C/L。
1.2.2 变分函数
变分计算是求泛函极值问题的一种数学方法。作为一种积分形式,如果将某一函数代入一个泛函中,则该泛函具有一个确定的数值。变分计算的主要问题是求函数f(x),使得函数的小变差δf(x)不至于改变原来的泛函。研究变差的计算并将其应用于有限元理论,会涉及线性代数、泛函分析和拓扑学原理等相关知识。
本节将介绍变分计算的基本理论,并由此介绍如何将泛函变分用于构造有限元模型。需要注意的是,上面各个方程的变分函数的用法,与应变能和最小势能原理在弹性理论和结构理论中的用法是相似的。
除含有一阶导数项的方程以外,上一节中其他控制方程的变分函数可以写成统一的形式。在本书的后续内容中,由于有限单元上的面积和材料弹性系数等被视作常数,因此上述方程中的相应项也将被视作常数。记为f=f(x),有
含有一阶导数项的方程则不一定有相应的变分函数。然而,为了得到上述方程的有限元模型,可以采用伪变分函数或拟变分函数来表示相应的控制微分方程。如果记式(1-16)中的C=C(x),则与其相应的拟变分函数为
对于本节所论述的控制微分方程,式(1-23)经过适当的变形即可得到相应方程的特征变分函数,而采用类似的分析过程却不能将式(1-24)经过适当变形导出式(1-16)的特征变分函数。但是,可以由式(1-24)导出式(1-16)的有限元模型。
在一般情况下,变分格式中包含边界条件。因而从这个角度来讲,式(1-23)的表示方式并不完善。然而,变分函数给出了函数与微分方程之间的对应关系,并成为初步研究有限元法的必要条件。
1.2.3 插值函数
有限元法的基本概念是连续函数可以近似地表示为离散模型。这里的离散模型是由一个或多个插值多项式组成的,而连续函数被分为有限段(或有限片、有限块),即有限个单元且每个单元由一个插值函数所定义,用以刻画单元在端点之间的状态。有限单元的端点称为节点。
1.2.4 形函数
形函数(在某些文献中也称为形状函数)常用字母N来表示,它通常为插值多项式的系数。在某个有限单元中,不同的节点有其各自的形函数,例如,某形函数在该节点的函数值为1,而在该单元上的其他节点处的函数值为0。插值多项式和形函数这两个概念经常交替使用。
1.2.5 刚度矩阵
刚度矩阵这个名词来源于结构分析,有限元法的最早期应用类似于矩阵的结构分析,用来描述力和位移之间的矩阵关系。目前,在提到刚度矩阵时,不再考虑其应用。温度与热通量之间的矩阵关系也称为刚度矩阵。
有限元法定义了两个刚度矩阵:单元刚度矩阵,对应于独立的某个单元;整体刚度矩阵,由所有单元刚度矩阵组合而成,定义的是整个系统的刚度矩阵。
1.2.6 连通性
连通性是指有限元模型汇总的一个单元与相邻单元的连接特征。在本章所论述的内容中,对于每个节点处有一个未知量的一维线性两节点单元而言,局部单元上的微分运算是主要的。上述单元的左端点的编号为1,右端点的编号为2。显然,整体有限元模型中的所有点不能都记为点1或点2。整体模型和局部模型之间通过连通度矩阵进行联系。如果整体模型中含有Ne1个有限单元,每个单元中含有Nnode个节点,则连通度矩阵的维数为Ne1×Nnode。如图1-1所示,整体模型中含有5个使用罗马字母表示的单元,而局部模型和整体模型之间通过一个如表1-1所示的5×2阶连通度矩阵进行联系。
图1-1 整体及局部有限元模型
表1-1 连通度矩阵
1.2.7 边界条件
在1.2.1节中,边界条件分为本质边界条件和自然边界条件两种,采用解析法可以得到类似于本章所介绍的二阶方程的解,但是这里需要计算两个积分常数。为了计算上述积分值,必须给出两个边界条件,这些边界条件通常在问题的一维定义域的两端分别给出。边界条件一般按照未知量的具体数学形式进行分类。
在数学术语中,本质边界条件也称为狄利克雷(Dirichlet)边界条件。由式(1-10),长度为L的等截面均匀直杆的一维定常状态的热传导问题,可以表示为如下的狄利克雷问题:
式中的两个边界条件都是针对温度给定的。
实际问题的诺依曼(Neumann)边界条件在两个边界点上都是给出的一阶导数值,这类问题称为诺依曼问题。对于热传导问题,其诺依曼边界条件给出的是通量所满足的条件,例如,联合式(1-25)和式(1-26)可得
这类边界条件给理论求解和数值求解都带来了不少困难,只有在给定某一温度的条件下,上述问题才是唯一可解的。本章不讨论类似于式(1-27)的边界条件的有限元分析。
第三类边界条件称为混合边界条件。这类边界条件相当于式(1-26)和式(1-27)的组合形式,是最常用的一类边界条件。实际上有两类混合边界条件。第一类混合边界条件是一个边界条件为本质边界条件,另一个边界条件为自然边界条件。第二类混合边界条件,如热传导方程的边界条件为
式中,h为对流换热系数,T∞为边界面以外介质的温度。上述边界条件表明边界通量与边界温度结合等于某一已知温度,第二类边界条件可以为本质边界条件和自然边界条件,也可以为类似于式(1-28)的形式。
1.2.8 圆柱坐标系中的问题
类似于热传导和静电分布等一维轴对称微分方程与式(1-10)和式(1-22)是相似的,与式(1-22)相对应的轴对称柱面坐标系下的电势方程为
由于此处的定义域对应于柱面边界的轴线,故面积为常数。式(1-29)可以简写为
与之相对应的变分函数为
式中,2πrdr代替了dV。
1.2.9 直接方法
直接方法通常用来叙述从矩阵结构分析到有限元常用概念的发展过程,其主要目的就是利用从材料力学推导出的直杆或梁单元问题的解法。例如,在材料力学中,在定常外力P的作用下,直杆上的轴向应力为σ=P/A,将以上定义与式(1-31)相结合,可以推导出直杆沿着轴向外力作用下的变形:
式中,μ是直杆的两端沿着轴向受到均匀外力P的作用下,长度为L的直杆的整体变形,由式(1-32)可以给出刚度矩阵。基于式(1-32)的推导过程一般被包括在桁架的应力分析中,这里每个桁架都可以被视作一个单元。外力P通常为桁架节点处的载荷。