張勝

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

1 引言

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

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

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

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

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

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

2 薄殼方程

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

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

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

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

2.1 Naghdi薄殼模型

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

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

2.2 殼厚度趨于零時(shí)殼模型解的漸近估計(jì)

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

3 有限元模型

4 有限元法的誤差分析

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

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