具有一致精度的薄殼方程的線性有限元方法

張勝

(韋恩州立大學，底特律 密歇根 48202)

1 引言

在Naghdi薄殼模型中，總應變能是彎曲應變能，切向延壓應變能和橫向剪切應變能的總和.在彎曲為主的變形中，殼的延壓剪切應變相對很小，當彎曲的殼無限變薄時延壓剪切應變趨向于零.如果把殼的變形限制在分片多項式構成的有限元函數(shù)空間里，當延壓剪切應變?yōu)榱銜r，有限元函數(shù)能表示的變形縮減為零.對很薄的殼而言，有限元解給出的變形遠小于實際變形.此所謂延壓剪切數(shù)值閉鎖.這個問題的根源在于分片多項式不能準確表達曲面的等距純彎曲變形.這種失敗的數(shù)值計算會導致工程師做出錯誤的關于殼結構強度的判斷.

數(shù)值閉鎖會發(fā)生在許多依賴參數(shù)的數(shù)學物理方程的科學計算中，這包括Timoshenko梁彎曲方程，Timoshenko-Naghdi拱變形方程，Reissner-Mindlin板彎曲方程，忽略了橫向剪切應變的 Koiter[1]薄殼方程，更一般的 Naghdi[2]殼方程和其他薄或細的構件的形變力學問題的模型.這些問題中的參數(shù)是構件的相對厚度.切向延壓閉鎖是Koiter殼模型數(shù)值計算中的一個主要問題，而橫向剪切閉鎖一直是Reissner-Mindlin板彎曲模型研究的中心問題.Naghdi薄殼模型則涉及這兩種閉鎖，而且殼的曲率使這兩種閉鎖耦合在一起無法分離.平板是曲率為零的特殊的殼，在此情況下Naghdi薄殼模型解耦成Reissner-Mindlin板方程和一個平板切向延壓的平面應力方程，對這兩者都有幾個成功的算法.但兩者簡單的結合無法產生有效的殼有限元.

對Timoshenko梁彎曲方程而言，在形成剛度矩陣的過程中只須用一個低精度的數(shù)值積分方法來計算剪切應變能便可消除剪切閉鎖，從而得到最優(yōu)階的一致精確的精度不隨梁的厚度變化的有限元方法.這個技術早已為結構工程師所知，其數(shù)學理論則要用混合有限元方法來建立[3].對于Reissner-Mindlin板，最成功的方法是基于對其變量的重組，把板方程解構成一個攝動過的Stokes方程和Poisson方程，組合其已有的有限元而得到的，參閱文獻[4]及其中的參考文獻.由于其特殊性，這些方法和理論無法用于薄殼模型.對于薄殼問題，自從有限元創(chuàng)立以來，工程力學界和數(shù)學界一直在不斷地努力，盡管有大量數(shù)值計算的工程文獻和很大的進展，可消除閉鎖的方法的數(shù)學理論還遠不如人意[5-8].幾個商用軟件都有各自的算法，但沒有一種方法是有數(shù)學基礎的或完全可靠的，有時甚至是失敗的.

間斷有限元(Discontinuous Galerkin簡稱DG[9])近年來得到了長足的發(fā)展.它給選擇有限元空間和自由度提供了更靈活的方法，在有些計算問題中，產生了高精度高效率算法.在理論上它有可能把五花八門的有限元納入一個統(tǒng)一的框架.不少人相信DG具有解決薄殼計算中閉鎖問題的潛力[10-11].文獻[12]分析了DG方法在解決Koiter殼計算中切向延壓閉鎖的問題優(yōu)勢.本文討論Naghdi薄殼模型的最低階混合DG方法，所用有限元函數(shù)均為分片線性函數(shù).對殼中面位移和法向纖維轉角用間斷函數(shù)(在與殼的自由邊界相臨的單元上須增加一些二次函數(shù))，而對輔助性的延壓應力張量和剪切應力向量用連續(xù)函數(shù).用Nitsche方法處理固支邊界，并繞過了混合有限元方法[13]常用的Babuˇska-Brezzi條件.如果用常規(guī)的線性有限元計算Naghdi方程，在以彎曲為主的殼變形問題中會有非常嚴重的閉鎖現(xiàn)象，致使數(shù)值結果完全無用.

這里的分析主要針對的是彎曲為主的殼變形問題.需要說明的是有些情況下，殼的變形是以切向延壓為主的，這時薄殼具有極高的承載能力.更多的殼變形是介于彎曲為主和延壓為主的中間形態(tài)，其承載能力高于彎曲的殼，但不如延壓殼抗載.殼變形屬于何種形態(tài)取決于薄殼曲面的形狀，加載方式和邊界支撐方式.例如，如果殼的中面是直紋面，部分邊界是一條直紋線，沿其固支，在橫向載荷作用下它的變形便是以彎曲為主.如果殼的中面是橢圓形的，沿整個邊界固支，不管如何加載，它的變形都以切向延壓為主.如果殼的中面是橢圓形的，沿部分邊界固支，其余部分自由，它的變形處于中間形態(tài).能夠避免數(shù)值閉鎖，從而對殼彎曲問題有效的算法是否適用于其它種類的殼變形是計算工程力學中的未曾解決的重大問題.本文亦無意做此嘗試.

本文結構如下，在第2節(jié)中，引進Naghdi殼模型，引入橫向剪切應力向量和切向延壓應力張量作為新變量，把殼方程寫成混合形式，并給出一些必要解的先驗估計.在第3節(jié)中，引入有限元模型.在第4節(jié)做誤差分析.在文中，C代表常數(shù)，其值可依賴于殼的曲率和其它幾何系數(shù)，殼的材料的Lamé系數(shù)，和有限元單元形狀的規(guī)則性有關，但與有限單元的尺寸和殼的厚度無關.用A?B來表示A≤CB.如果A?B和A?B都成立，寫成A?B.用上標表示向量和張量的反變分量，下標表示協(xié)變分量.除?外，希臘字母上下標在{1，2}中取值.拉丁字母在{1，2，3}中取值.也采用關于重復上下標的Einstein加法規(guī)則，和Sobolev空間中的常用記號.具有協(xié)變分量uα或反變分量ξα的向量將分別由粗體字母u或ξ表示.具有分量Mαβ的張量將簡稱為M.

2 薄殼方程

協(xié)變微分的乘法規(guī)則，如 (σαλuλ)|β=σαλ|βuλ+σαλuλ|β也是成立的.

映射φ是Ω和之間的一對一對應關系，它把子域τ?Ω映射到子區(qū)域

曲線段e?映射成曲線段=φ(e).在殼中面上定義的函數(shù)f將通過映射φ與Ω上定義的函數(shù)認同，并用相同的符號表示.因此f(φ(xα))=f(xα).若無進一步解釋，波浪號表示曲面上的量或運算，沒有波浪號則表示在平面域Ω上操作.需要使用曲面上的格林公式，反復進行分部積分.對曲面子域，用表示與曲面相切的邊界?=φ(?τ)的單位法向量.設nαeα是R2中?τ的單位法向量.這里的eα是R2中的基向量.對于向量場fα，格林公式如下:

在這里和下文中，為簡單計，忽略了積分中的微分元素.第一個積分是關于曲面的面積元的，表達成平面區(qū)域τ上的常規(guī)積分則是∫第二個積分是根據?的弧長取的.最后一個是根據弧長?τ.將使用這樣一個事實，在?τ的直線部分上，nα是常數(shù)，而α通常是沿?變化的.

2.1 Naghdi薄殼模型

Naghdi殼模型[2]使用殼中面的切向位移u=uαaα，法向位移wa3和法向纖維旋轉θ=θαaα作為主要變量.用這樣一組主要變量，彎曲應變，切向延壓應變和橫向剪切應變可表示如下:

這個混合模型是有限元方法的基礎，它的解由十個定義在二維區(qū)域Ω上的函數(shù)組成.

2.2 殼厚度趨于零時殼模型解的漸近估計

引用文獻[17-18]中的兩個結果，得到一些關于Naghdi殼模型的解隨殼厚度變化的漸近行為的有用估計，以分析有限元模型.本小節(jié)中的符號獨立于本文的其余部分.Naghdi殼模型(2.7)可寫成如下算子方程(2.12).假設H，U，V是Hilbert空間，A和B是分別從H到U和V的線性連續(xù)算子.假設

3 有限元模型

4 有限元法的誤差分析

有限元模型(3.9)的穩(wěn)定性和它與Naghdi殼模型(2.10)的相容性保證了有限元解的最佳逼近性.有限元法的誤差分析簡化成了一個逼近論問題.對Naghdi殼方程的解，可構造有限元空間的插值函數(shù)，從而證明如下定理，它是本文的主要結果.這個定理的證明和相關數(shù)值驗證可參考文獻[12]和它引用的文獻.

這個極限是個非零常數(shù).定理4.1意味著有限元解的相對誤差具有最優(yōu)階的精度.本方法解決閉鎖問題的效率從定理中不等式右側第一個括號中的項可以看出.如果用常規(guī)有限元方法這括號中的項將變成[1+?-1]，當?→0這個系數(shù)會被無限放大，使計算結果對很薄的殼失效.如果殼變形不是彎曲主導的，則有‖(θ0，u0，w0)‖Hh=0，定理 4.1并不意味著有限元解在相對誤差下具有任何精度.