張琪昌，趙 彬，王 煒

(1. 天津大學(xué)機械工程學(xué)院，天津 300072；2. 天津大學(xué)天津市非線性動力學(xué)與混沌控制重點實驗室，天津 300072)

非對稱截面環(huán)狀Kirchhoff彈性細(xì)桿的拓?fù)錁?gòu)型分析

張琪昌1,2，趙 彬1,2，王 煒1,2

(1. 天津大學(xué)機械工程學(xué)院，天津 300072；2. 天津大學(xué)天津市非線性動力學(xué)與混沌控制重點實驗室，天津 300072)

為了揭示非對稱截面環(huán)狀Kirchhoff彈性細(xì)桿受到外力(矩)作用時其拓?fù)錉顟B(tài)的變化情況，引入了等效初始扭矩的概念，利用勢能密度函數(shù)與哈密頓函數(shù)形式上的相似性，建立了其與系統(tǒng)哈密頓函數(shù)、初始自連接數(shù)之間的定量關(guān)系．以環(huán)狀DNA分子為背景，對扭矩模型進(jìn)行了非線性動力學(xué)分析，并利用解析方法中的待定固有頻率法，分析了主要參數(shù)對扭矩模型穩(wěn)定構(gòu)型的影響．最后，利用等效剛度得到相應(yīng)狀態(tài)下曲率的表達(dá)式，為進(jìn)一步認(rèn)識和描述DNA環(huán)受到生物酶作用下的穩(wěn)定構(gòu)型提供了新的理論途徑．

Kirchhoff；彈性細(xì)桿；環(huán)狀DNA分子；待定固有頻率法；非線性動力學(xué)；等效初始扭矩

關(guān)于環(huán)狀彈性細(xì)桿的研究，目前所采用的模型是在Kirchhoff假設(shè)的基礎(chǔ)上，將直桿兩端的邊界條件進(jìn)行等效處理，從而衍生出環(huán)狀形態(tài)的桿[1-4]．Shi等[4]在彈性桿Kirchhoff方程基礎(chǔ)上，以彈性桿的曲率和撓率為基本變量，建立了圓截面彈性桿的Schr?dinger模型，并獲得了其橢圓函數(shù)形式的近似解析解．劉延柱等[5-7]基于彈性桿非線性力學(xué)的經(jīng)典理論，研究了以分子生物學(xué)為背景的彈性桿平衡和穩(wěn)定性問題，發(fā)展了以歐拉角為坐標(biāo)的參數(shù)理論體系，在靜、動力學(xué)方面都獲得了大量很有啟示意義的結(jié)論．文獻(xiàn)[8-9]中引入復(fù)剛度概念，推導(dǎo)出了非圓截面彈性桿的Schr?dinger方程，并同樣獲得了其橢圓函數(shù)形式的近似解析解．王煒等[10]在考慮DNA彈性細(xì)桿截面非對稱性特征的基礎(chǔ)上，構(gòu)造出復(fù)數(shù)形式的Kirchhoff系統(tǒng)，利用復(fù)變量力矩設(shè)解形式，并對Kirchhoff高維系統(tǒng)進(jìn)行有效簡化，得到了以外力矩軸向分量為單獨變量的二階常微分方程．

眾所周知，DNA是由兩條超長堿基鏈按一定順序配對排列而成，由于細(xì)胞的空間約束性，現(xiàn)實中的DNA不存在直桿狀態(tài)，DNA鏈必須在特定的外力及外力扭矩作用下發(fā)生纏繞、彎曲，形成超螺旋狀態(tài)，基于此研究彈性桿拓?fù)錁?gòu)型及其發(fā)展演變，將有助于揭示其生物學(xué)本質(zhì)．文獻(xiàn)[1-2]中提出了一種基于能量最小原理與解析方法的判定方法，該方法可以判定環(huán)狀桿構(gòu)型的分岔現(xiàn)象．與數(shù)值方法中通常采用的Monta Carlo理論相比，在所得結(jié)果基本一致的前提下，定量分析不僅計算耗時短，而且還可以有效地反映拓?fù)鋮?shù)之間的相互關(guān)系．正是基于上述考慮，作者提出了用連接數(shù)和纏繞數(shù)來描述彈性桿構(gòu)型的新方法．

為深入分析外力(矩)作用對彈性桿構(gòu)型的影響，本文首先引入了等效初始扭矩的概念，運用解析方法，充分分析了其對系統(tǒng)參數(shù)和扭矩模型參數(shù)的影響．再對文獻(xiàn)[10]所提出的扭矩模型進(jìn)行進(jìn)一步探討，結(jié)合非線性動力學(xué)待定固有頻率法[11-12]，求出了形式更為直觀的近似解析解，利用等效剛度的概念，求出了彈性桿撓性線曲率的表達(dá)式．結(jié)合環(huán)狀桿的邊界條件，分析了不同力矩作用下桿的姿態(tài)特性變化，畫出了相應(yīng)狀態(tài)下?lián)闲跃€的實際效果圖．

1 等效初始扭矩

任意一個初始扭率為0的環(huán)狀彈性桿姿態(tài)均可由這樣的形式得到:如圖1所示，在一個撓性線初始扭率為0的非圓截面彈性桿上，設(shè)定點O為弧長原點(弧長s變化范圍為0～L，L為桿長)，在s=0和s=L處兩截面分別施以大小相等、方向相反的等效初始扭矩對M30(沿彈性桿中心線切線方向)和主矢對F，使得桿在兩端界面發(fā)生扭轉(zhuǎn)的同時，保持兩界面的距離始終為無限小，當(dāng)達(dá)到既定的狀態(tài)時，將兩端重新粘結(jié)在一起，再次形成閉合的回環(huán)．非對稱截面桿的截面如圖2所示，其中Rxy和Rz分別表示垂直軸線上的長短半軸，且Rxy∶Rz≈1．

由于在這個過程中，主矢F的作用只限于控制s=0和s=L兩截面之間的距離，使二者不至于分離太遠(yuǎn)，而扭矩M30才對該環(huán)狀彈性桿輸入有效能量．為保證首尾截面上同一撓性線的首尾端曲率以及曲率變化情況一致，初始扭矩M30的做功扭角必須滿足φ=2nπ，n為正整數(shù)．在上述過程中，主矢F只是起到端面軸向固定作用，并沒有能量輸入，大小可以忽略．而初始扭矩M30在做功完成之后使彈性桿的首尾截面相對扭角為φ，桿自身作為傳導(dǎo)兩端力矩作用的載體，其內(nèi)部的應(yīng)力大小直接由M30決定．

圖1 環(huán)狀桿受到等效初始扭矩作用Fig.1 Effect of equivalent initial torque on ring rod

圖2 非對稱截面桿的截面(Rxy∶Rz≈1)Fig.2 Cross section of asymmetric elastic rod(Rxy∶Rz≈1)

綜上所述，由于初始扭矩方向為沿截面法線方向，不具有其他方向的分量，且首尾截面沒有發(fā)生相對位移和扭角，主矢F(包括其軸向分量F3)對系統(tǒng)的做功為0，可以將哈密頓函數(shù)H的定義式簡化為

式中:A、B均為彈性桿抗彎剛度；C為抗扭剛度．由式(1)可以看出，此時等效初始扭矩M30的大小，決定了彈性桿再次閉合后封閉系統(tǒng)中的能量，據(jù)此，可以由M30對封閉環(huán)狀桿的狀態(tài)進(jìn)行標(biāo)定．

根據(jù)勢能最小原理[5]，彈性桿的實際平衡狀態(tài)對應(yīng)于總彈性勢能的極小值，即總勢能為0的狀態(tài)，故有

式中:δE為總勢能的一階變分；Γ為勢能密度函數(shù)．由于在忽略主矢F的前提下，勢能密度函數(shù)Γ=H ，首尾端面相對虛角位移為δφ= φ．式(2)中第2個等號兩邊對弧長s求導(dǎo)，再考慮自連接數(shù)的定義，可知初始自連接數(shù)Lk0與M30的關(guān)系為

生物學(xué)上，DNA的扭轉(zhuǎn)、剪切以及聚合等過程是由各種各樣的酶來主導(dǎo)并實施的，例如DNA限制性內(nèi)切酶、解旋酶、連接酶等，直接定量研究這些酶對DNA作用難度較大，而等效初始扭矩則可以在標(biāo)定外界酶對DNA環(huán)的能量輸入、預(yù)測DNA環(huán)在受一定扭矩后的穩(wěn)定構(gòu)型以及確定環(huán)狀桿在某一狀態(tài)下的數(shù)學(xué)模型等方面發(fā)揮重要作用．這也是本文研究環(huán)狀彈性桿的拓?fù)鋮?shù)進(jìn)而開展動力學(xué)分析的前提條件．

2 動力學(xué)分析

非對稱截面Kirchhoff桿的模型[10]為

式中:M3為扭矩沿弧長s的分布函數(shù)；M′和M′′分別為M3對弧長s的一階和二階導(dǎo)數(shù)．

考慮實際參數(shù)意義可知，抗彎、抗扭剛度以及哈密頓函數(shù)值均為正值，所以式(4)中，等號右邊關(guān)于扭矩M3一次方項的系數(shù)滿足

該系數(shù)的相反數(shù)對應(yīng)于常見振動微分方程中回復(fù)力的系數(shù)，即固有頻率的平方，所以，該系統(tǒng)的平衡點可以在圖3中得到定性的描述．利用穩(wěn)定性理論可以判斷出中間的坐標(biāo)原點為鞍點，周圍的軌線是不穩(wěn)定的，而原點左右兩邊各分布一個中心點(R-，0)和(R+，0)，其周圍分布有穩(wěn)定的軌線環(huán)1和環(huán)2，說明這兩點附近的解是穩(wěn)定的，而相應(yīng)扭矩M3的周期變化幅值若超過一定范圍，則表示其軌線已經(jīng)達(dá)到甚至越過原點，形成跨度更大的環(huán)3和環(huán)4，扭矩M3出現(xiàn)零值，甚至跨越遠(yuǎn)點，出現(xiàn)負(fù)值，很顯然這會使桿的姿態(tài)變得不穩(wěn)定．

圖3 扭矩模型的平衡點及相應(yīng)的相圖Fig.3 Equilibrium points and corresponding phase diagrams of torque model

根據(jù)實際意義，只有扭矩M3在正值區(qū)域內(nèi)做周期運動時才是符合工程實際情況的，所以本文選擇對右平衡點進(jìn)行分析．由原振動方程可以知道，該平衡點(R+，0)的坐標(biāo)其中

為便于計算和書寫，設(shè)

其中

利用待定固有頻率法求得該振動系統(tǒng)的近似解析解為

式中:a為振幅；ω1為待定圓頻率，二者關(guān)系為

幅值a隨待定圓頻率ω1的變化曲線如圖4所示，可見a隨ω1的增大而增大．環(huán)狀桿的邊界條件為

式中κ(s)為曲率沿弧長的分布函數(shù)．結(jié)合式(17)和式(3)，可以得出

可見當(dāng)桿長以及桿的材料特性不變時，環(huán)狀彈性桿扭矩變化的圓頻率與初始等效扭矩大小成正比關(guān)系．

綜合式(13)～式(17)，可以確定出扭矩的近似解析解，這意味著在運用解析方法的基礎(chǔ)上，彈性桿穩(wěn)定構(gòu)型受各個參數(shù)變化的影響而發(fā)生的變化情況可以動態(tài)的形式得出．同時，由于該扭矩模型可進(jìn)一步發(fā)展和完善，研究者根據(jù)模型參數(shù)情況，可以更加準(zhǔn)確、清晰地判斷其對桿構(gòu)型的影響情況．

利用非對稱截面等效剛度Beff的定義[10]可知

根據(jù)本構(gòu)關(guān)系

可以得出

圖4 幅值隨待定圓頻率的變化曲線Fig.4 The ω1-a curve

從而得到了非對稱截面Kirchhoff桿撓性線曲率的表達(dá)式．上述討論的右平衡點附近扭矩的周期變化區(qū)域始終為正，這也恰恰符合曲率平方的取值范圍.

3 數(shù)值模擬

本節(jié)將以第2節(jié)得到的解析解為基礎(chǔ)，帶入具體數(shù)值，分析等效初始扭矩對環(huán)狀桿穩(wěn)定構(gòu)型的影響．

對彈性桿模型的物理參數(shù)進(jìn)行取值:抗彎剛度分別為A=0.13,erg·cm和B=0.12,erg·cm，抗扭剛度C=0.10,erg·cm，桿長L=20,cm，Rxy∶Rz=1.3∶1.2．將其代入式(1)～式(3)以及式(7)～式(17)，并將坐標(biāo)原點移至原始系統(tǒng)的原點處，即M3=3M++ R. 取M30=1.6π N·cm，相應(yīng)自連接數(shù)0kL=4，可得

再將式(22)代入式(21)，就可以進(jìn)一步得到曲率平方的表達(dá)式為

取M30=2π N·cm，相應(yīng)自連接數(shù)0kL=5，可以得到扭率近似解以及相應(yīng)的撓性線曲率表達(dá)式為

利用Mathematica軟件，分別畫出M30=1.6π N·cm(0kL=4)和M30=2π N·cm(0kL=5)時扭矩變化相圖和相應(yīng)的隨弧長變化的歷程圖，如圖5和圖6所示．可以看出，等效初始扭矩越大，形成穩(wěn)定狀態(tài)時它在環(huán)狀桿上功也越大，扭矩沿弧長的分布就會越來越密集，幅值也越來越大；從圖6所示的相圖上看，穩(wěn)定環(huán)的半徑也越來越大．

由此可見，DNA的雙螺旋結(jié)構(gòu)在復(fù)制完成時，在酶的扭矩作用下，堿基鏈所代表的撓性線曲率振動幅值不斷變大，封閉環(huán)所蘊含的能量也越來越大．圖3所示相圖中，圓環(huán)與橫軸的左交點也逐漸靠近原點．由于(0，0)點為鞍點，所以撓性線曲率的最小值越接近0值，曲桿狀態(tài)中出現(xiàn)直線狀態(tài)越不穩(wěn)定，當(dāng)堿基鏈的曲率在受到外界干擾時極易跳躍到負(fù)值，即正反向扭轉(zhuǎn)同時出現(xiàn)在一個DNA環(huán)中，桿的狀態(tài)即出現(xiàn)失穩(wěn)．

圖5 扭矩周期變化的相圖Fig.5 Thecurve

圖6 扭矩隨弧長變化的歷程Fig.6 The s-M3curve

同樣，根據(jù)式(23)和(25)所示的曲率變化情況，畫出環(huán)狀桿撓性線效果圖，見圖7和圖8．

圖7 M30=1.6π N·cm對應(yīng)的撓性線效果圖Fig.7 Flexibility of the rod with M30=1.6π N·cm

圖8 M30=2π N·cm對應(yīng)的撓性線效果圖Fig.8 Flexibility of the rod with M30=2π N·cm

4 結(jié) 論

(1) 根據(jù)勢能密度函數(shù)與哈密頓函數(shù)的形式相似性，引入等效初始扭矩的概念．結(jié)果表明，環(huán)狀彈性桿作為一個封閉系統(tǒng)，其哈密頓函數(shù)的數(shù)值大小與等效初始扭矩的平方成正比關(guān)系，而初始自連接數(shù)與其一次方成正比關(guān)系，利用該參數(shù)，可以確定在DNA環(huán)復(fù)制、纏繞、螺旋和解螺旋等過程中，生物酶對封閉環(huán)的能量輸入作用，以及對穩(wěn)態(tài)構(gòu)型形成機制的影響等功能．

(2) 運用非線性動力學(xué)思想對彈性桿扭矩模型進(jìn)行了分析，系統(tǒng)分析了非對稱截面Kirchhoff彈性桿扭矩沿模型所有平衡點的情況，并重點考慮和分析了符合實際情況的中心平衡點，以解析法為基礎(chǔ)，分析了等效初始扭矩對該動力學(xué)模型的影響．在中心平衡點附近，當(dāng)桿長以及桿的材料特性不變時，等效初始扭矩越大，環(huán)狀彈性桿扭矩變化的圓頻率越大，扭矩和撓性線曲率分布的變化范圍和幅值也就越大，在相圖中的運動軌線就越接近原點，正反繞向撓性線的同時存在，將導(dǎo)致DNA雙堿基鏈在扭轉(zhuǎn)過程中極易出現(xiàn)失穩(wěn)情況．

(3) 本文分析過程及結(jié)論體現(xiàn)了解析方法的優(yōu)勢，即以直觀的方式表現(xiàn)出參數(shù)變化對彈性桿構(gòu)型的動態(tài)、連續(xù)影響．相圖中平衡點穩(wěn)定的中心特性也使得方程解析解的穩(wěn)定性更具說服力，當(dāng)然這也體現(xiàn)了李雅普諾夫理論無需求解即可判定解的穩(wěn)定性的優(yōu)點．相比數(shù)值方法來說，該解析表達(dá)式的形式更為直觀簡便，計算規(guī)模卻要小很多．本文結(jié)果同樣與Monta Carlo積分得出的結(jié)果相一致，符合能量最小原理的要求．本文分析過程不僅對于非對稱截面彈性桿的研究適用，對其他以Kirchhoff方程為基本模型的彈性桿研究也同樣具有重要的應(yīng)用價值．

［1］ Hoffman K A，Manning R S，Maddocks J H. Link，twist，energy，and the stability of DNA minicircles[J]. Biopolymers，2003，70(2):145-157.

［2］ Furrer P B，Manning R S，Maddocks J H. DNA rings with multiple energy minima[J]. Biophysical Journal，2000，79(1):116-136.

［3］ Stump D M，Watson P J，F(xiàn)raser W B. Mathematical modelling of interwound DNA supercoils[J]. Journal ofBiomechanics，2000，33(4):407-413.

［4］ Shi Yaoming，Hearst J E. The Kirchhoff elastic rod，the nonlinear Schr?dinger equation，and DNA supercoiling[J]. Journal of Chemical Physics，1994，101(6): 5186-5200.

［5］ 劉延柱. 彈性細(xì)桿的非線性力學(xué)——DNA力學(xué)模型的理論基礎(chǔ)[M]. 北京:清華大學(xué)出版社，2006.

Liu Yanzhu. Nonlinear Mechanics——Thin Elastic Rod-Theoretical Basis of Mechanical Model of DNA[M]. Beijing:Tsinghua University Press，2006(in Chinese).

［6］ Liu Yanzhu，Xue Yun. Stability analysis of helical rod based on exact Cosserat model[J]. Applied Mathematics Mechanics，2011，32(5):603-612.

［7］ Liu Yanzhu，Sheng Liwei. Stability and vibration of a helical rod constrained by a cylinder [J]. Acta Mechanica Sinica，2007，23(2):215-219.

［8］ 薛 紜，劉延柱. Kirchhoff彈性直桿在力螺旋作用下的穩(wěn)定性[J]. 物理學(xué)報，2009，58(10):6437-6442.

Xue Yun，Liu Yanzhu. Stability of a straight Kirchhoff elastic rod under the force screws[J]. Acta Physica Sinica，2009，58(10):6437-6442(in Chinese).

［9］ Xue Yun，Liu Yanzhu，Chen Liqun. Schr?dinger’s equation for the Kirchhoff elastic rod with non-circular cross section [J]. Chinese Physics，2004，13(6):794-797.

［10］ 王 煒，張琪昌，靳 剛. 非對稱截面Kirchhoff彈性細(xì)桿模型簡化方法研究[J]. 物理學(xué)報，2011，61(6):0646021-0646026.

Wang Wei，Zhang Qichang，Jin Gang. The analytical reduction of the Kirchhoff thin elastic rod model with asymmetric cross section[J]. Acta Physica Sinica，2011，61(6):0646021-0646026(in Chinese).

［11］ Leung A Y T，Zhang Q C. Complex normal form for strongly non-linear vibration systems exemplified by Duffing-van der Pol equation[J]. Journal of Sound and Vibration，1998，213(5):907-914.

［12］ Wang Wei，Zhang Qichang，F(xiàn)eng Jingjing. Global bifurcations of strongly nonlinear oscillator induced by parametric and external excitation[J]. Science China，2011，54(8):1986-1991.

Topological Configuration of Kirchhoff Thin Elastic Ring Rod with Asymmetric Cross Section

Zhang Qichang1,2，Zhao Bin1,2，Wang Wei1,2
(1. School of Mechanical Engineering，Tianjin University，Tianjin 300072，China；2. Tianjin Key Laboratory of Nonlinear Dynamics and Chaos Control，Tianjin University，Tianjin 300072，China)

Considering the forces and torques acted on the Kirchhoff thin elastic ring rod with asymmetric cross section，equivalent initial torque was introduced as a new concept to reveal the quantitative relationship between the state properties and parameters of torque dynamic model，which determines the configuration of the ring rod. As there exists formally a similarity between the potential energy density function and the Hamiltonian function，the quantitative relationship was obtained，which connects the equivalent initial torque，Hamiltonian function and initial link number. With DNA ring molecule as the physical background，this paper made the nonlinear dynamic analysis of torque model，and used undecided fundamental frequency method to obtain the asymptotic expressions for steady state periodic solutions around the practical equilibrium point. When the equivalent initial torque took different values，the graph of arc length-torque and phase portrait were drawn. The curvature could be expressed by the torque function with the equivalent bending stiffness considered，which provides a new methodological way to understand and describe the stable topological configuration of DNA ring under the influences of biological enzymes.

Kirchhoff；thin elastic rod；DNA ring molecule；undecided fundamental frequency method；nonlinear dynamic；equivalent initial torque

O322

A

0493-2137(2013)12-1089-06

DOI 10.11784/tdxb20131206

2012-05-11；

2012-09-02.

國家自然科學(xué)基金資助項目(11072168)；高等學(xué)校博士學(xué)科點專項科研基金資助項目(新教師類)(20100032120006).

張琪昌（1959— ），男，教授.

張琪昌，qzhang@tju.edu.cn.

時間：2013-11-08.

http://www.cnki.net/kcms/detail/12.1127.N.20131108.1610.017.html.