非對稱截面環狀Kirchhoff彈性細桿的拓撲構型分析

張琪昌，趙 彬，王 煒

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

非對稱截面環狀Kirchhoff彈性細桿的拓撲構型分析

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

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

為了揭示非對稱截面環狀Kirchhoff彈性細桿受到外力(矩)作用時其拓撲狀態的變化情況，引入了等效初始扭矩的概念，利用勢能密度函數與哈密頓函數形式上的相似性，建立了其與系統哈密頓函數、初始自連接數之間的定量關系．以環狀DNA分子為背景，對扭矩模型進行了非線性動力學分析，并利用解析方法中的待定固有頻率法，分析了主要參數對扭矩模型穩定構型的影響．最后，利用等效剛度得到相應狀態下曲率的表達式，為進一步認識和描述DNA環受到生物酶作用下的穩定構型提供了新的理論途徑．

Kirchhoff；彈性細桿；環狀DNA分子；待定固有頻率法；非線性動力學；等效初始扭矩

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

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

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

1 等效初始扭矩

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

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

圖1 環狀桿受到等效初始扭矩作用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(包括其軸向分量F3)對系統的做功為0，可以將哈密頓函數H的定義式簡化為

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

根據勢能最小原理[5]，彈性桿的實際平衡狀態對應于總彈性勢能的極小值，即總勢能為0的狀態，故有

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

生物學上，DNA的扭轉、剪切以及聚合等過程是由各種各樣的酶來主導并實施的，例如DNA限制性內切酶、解旋酶、連接酶等，直接定量研究這些酶對DNA作用難度較大，而等效初始扭矩則可以在標定外界酶對DNA環的能量輸入、預測DNA環在受一定扭矩后的穩定構型以及確定環狀桿在某一狀態下的數學模型等方面發揮重要作用．這也是本文研究環狀彈性桿的拓撲參數進而開展動力學分析的前提條件．

2 動力學分析

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

式中:M3為扭矩沿弧長s的分布函數；M′和M′′分別為M3對弧長s的一階和二階導數．

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

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

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

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

為便于計算和書寫，設

其中

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

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

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

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

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

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

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

根據本構關系

可以得出

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

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

3 數值模擬

本節將以第2節得到的解析解為基礎，帶入具體數值，分析等效初始扭矩對環狀桿穩定構型的影響．

對彈性桿模型的物理參數進行取值:抗彎剛度分別為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)，并將坐標原點移至原始系統的原點處，即M3=3M++ R. 取M30=1.6π N·cm，相應自連接數0kL=4，可得

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

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

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

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

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

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

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

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

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

4 結 論

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

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

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

［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，Fraser 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］ 劉延柱. 彈性細桿的非線性力學——DNA力學模型的理論基礎[M]. 北京:清華大學出版社，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彈性直桿在力螺旋作用下的穩定性[J]. 物理學報，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彈性細桿模型簡化方法研究[J]. 物理學報，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，Feng 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.

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

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

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

時間：2013-11-08.

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