基于祖沖之類方法的多體動力學方程保能量保約束積分

吳鋒　高強　鐘萬勰

摘要：針對一類多體動力學問題導出的微分代數方程，提出一種保能量、保約束的算法.該算法基于祖沖之類方法和歐拉中點保辛差分，利用祖沖之類方法保證在時間格點上精確滿足約束方程，避免約束違約問題；并進一步證明該算法在時間格點上可以精確保能量.數值算例進一步驗證該算法的可靠性.

關鍵詞：多體動力學方程； 微分代數方程； 保辛； 祖沖之

中圖分類號： O313.7

文獻標志碼： B

0引言

近年來，多體動力學研究已逐漸成為研究熱點.一方面是由于多體動力學建模得到的往往是強非線性方程，分析困難；另一方面則是因為多體動力學問題具有十分廣泛的運用背景[12]，如航空航天器、車輛和機器人等.目前，多體動力學的建模方法可以分為相對坐標法和絕對坐標法2類.當用絕對坐標法時，得到的往往是微分代數方程組.關于微分代數方程組的求解，目前已有一些進展.

許多學者認為，數值積分過程中約束方程的違約是造成積分困難的重要原因之一.原亮明等[3]把位移約束方程按照泰勒級數展開，與動力學方程組合進行迭代求解；趙維加等[4]用泰勒級數把約束方程展開，根據積分步長提出一種能對約束誤差自動修正的小擾動違約穩定方法；付士慧等[5]對具有完整、定常約束的多體系統，通過修改的帶乘子的拉格朗日正則形式方程，給出一種違約修正方法；戈新生等[6]提出一種基于完全笛卡爾坐標的多體系統微分代數方程符號線性化方法.在諸多研究中，文獻[7]中所提的方法能很好地解決約束方程違約問題.他們利用保辛的時間有限元結合祖沖之類方法的思想：保辛方法可以保證長時間積分的精度，祖沖之類方法可以解決約束方程違約這一問題.但是，當時的文獻中還沒有正式提出祖沖之類方法這一名詞，直到在最近出版的《經典力學辛講》[8]中，祖沖之類方法才被正式提出.祖沖之類方法指出，約束方程不必處處滿足，只要在時間格點處滿足即可.依據這一思想，使得微分代數方程的求解格式簡單，而且無約束違約問題.

對于非線性系統的差分格式要保辛，而國外論文提出，對于不可積系統，保辛和守恒難以同時達成[9]的問題，闡明保辛則能量不能守恒，能量守恒就不能保辛的兩難命題.實際上，文獻[10]通過引入含參變量的近似Hamilton系統，并以此為基礎利用保辛攝動的思想，提出一種Hamilton系統的保辛守恒積分算法，實現即保辛又保能量的算法.

圖7為能量的相對誤差，其中計算相對誤差所參考的真實能量H=34.868 J；圖8為約束最大相對誤差.由圖7和8可知，在積分1 000 s的時間內，本文算法計算結果的能量相對誤差和約束相對誤差都相當小，其中能量的相對誤差為10-13數量級，而約束誤差則為10-15數量級，說明本文算法完全是既保約束又保能量的，沒有約束違約問題，計算結果驗證理論分析結果.由于本算例雖然是在空間運動的雙擺，但是給定的初始位移和初始速度設定其實際運動只在一個平面內.圖9和10分別給出質點1的運動軌跡和質點2相對于質點1的運動軌跡.圖9和10中分別給出的3幅圖是從3個不同角度繪制的運動軌跡，其中，圖9（a）為xy平面內的運動軌跡，可見在計算1 000 s的時間區段內，質點1仍然精確地維持在一個平面內運動；同樣圖10（a）亦可驗證這一點.這進一步說明本文方法既保約束又保能量的優點，因為不守恒算法計算時，常常是在最初幾秒可以維持在同一平面；隨著計算繼續，誤差累計，運動軌跡往往躍出平面外.

3結束語

綜合運用文獻[78]提出的祖沖之類方法和文獻[10]提出的保辛守恒思想以及歐拉中點保辛差分格式，建立求解多體動力學方程的既保約束又保能量的算法.該方法在格點處可以嚴格滿足位移約束方程和能量守恒條件，而在非格點處，約束條件和能量方程可以近似滿足，避免約束違約的問題.數值算例表明，本文方法的計算結果滿意.

參考文獻：

[1]胡繼云. 建立多剛體系統動力學方程的坐標變換法及其應用[D]. 重慶： 重慶大學， 2004.

[2]田強. 基于絕對節點坐標方法的柔性多體系統動力學研究與應用[D]. 武漢： 華中科技大學， 2009.

[3]原亮明， 王成國， 劉金朝， 等. 一種求解多體系統微分代數方程的拉格朗日乘子方法[J]. 中國鐵道科學， 2001， 22（2）： 5154.

[4]趙維加， 潘振寬， 王藝兵. 多體系統動力學微分/代數方程約束誤差小擾動自我穩定方法[J]. 應用數學和力學， 2000， 21（1）： 9498.

[5]付士慧， 王琪. 多體系統動力學方程違約修正的數值計算方法[J]. 計算力學學報， 2007， 24（1）： 4449.

[6]戈新生， 趙維加， 陳立群. 基于完全笛卡爾坐標的多體系統微分代數方程符號線性化方法[J]. 工程力學， 2004， 21（4）： 106111.

[7]鐘萬勰， 高強. 約束動力系統的分析結構力學積分[J]. 動力學與控制學報， 2006， 4（3）： 193200.

[8]鐘萬勰， 高強， 彭海軍. 經典力學辛講[M]. 大連： 大連理工大學出版社， 2013.

[9]ZHONG G， MARSDEN J E. LiePoisson HamiltonJacobi theory and LiePoisson integrators[J]. Phys Lett A， 1988， 133（3）： 134139.

[10]高強， 鐘萬勰. Hamilton系統的保辛守恒積分算法[J]. 動力學與控制學報， 2009， 7（3）： 193199.

[11]邢譽峰， 楊蓉. 動力學平衡方程的Euler中點辛差分求解格式[J]. 力學學報， 2007， 39（1）： 100105.

（編輯武曉英）

摘要：針對一類多體動力學問題導出的微分代數方程，提出一種保能量、保約束的算法.該算法基于祖沖之類方法和歐拉中點保辛差分，利用祖沖之類方法保證在時間格點上精確滿足約束方程，避免約束違約問題；并進一步證明該算法在時間格點上可以精確保能量.數值算例進一步驗證該算法的可靠性.

關鍵詞：多體動力學方程； 微分代數方程； 保辛； 祖沖之

中圖分類號： O313.7

文獻標志碼： B

0引言

近年來，多體動力學研究已逐漸成為研究熱點.一方面是由于多體動力學建模得到的往往是強非線性方程，分析困難；另一方面則是因為多體動力學問題具有十分廣泛的運用背景[12]，如航空航天器、車輛和機器人等.目前，多體動力學的建模方法可以分為相對坐標法和絕對坐標法2類.當用絕對坐標法時，得到的往往是微分代數方程組.關于微分代數方程組的求解，目前已有一些進展.

許多學者認為，數值積分過程中約束方程的違約是造成積分困難的重要原因之一.原亮明等[3]把位移約束方程按照泰勒級數展開，與動力學方程組合進行迭代求解；趙維加等[4]用泰勒級數把約束方程展開，根據積分步長提出一種能對約束誤差自動修正的小擾動違約穩定方法；付士慧等[5]對具有完整、定常約束的多體系統，通過修改的帶乘子的拉格朗日正則形式方程，給出一種違約修正方法；戈新生等[6]提出一種基于完全笛卡爾坐標的多體系統微分代數方程符號線性化方法.在諸多研究中，文獻[7]中所提的方法能很好地解決約束方程違約問題.他們利用保辛的時間有限元結合祖沖之類方法的思想：保辛方法可以保證長時間積分的精度，祖沖之類方法可以解決約束方程違約這一問題.但是，當時的文獻中還沒有正式提出祖沖之類方法這一名詞，直到在最近出版的《經典力學辛講》[8]中，祖沖之類方法才被正式提出.祖沖之類方法指出，約束方程不必處處滿足，只要在時間格點處滿足即可.依據這一思想，使得微分代數方程的求解格式簡單，而且無約束違約問題.

對于非線性系統的差分格式要保辛，而國外論文提出，對于不可積系統，保辛和守恒難以同時達成[9]的問題，闡明保辛則能量不能守恒，能量守恒就不能保辛的兩難命題.實際上，文獻[10]通過引入含參變量的近似Hamilton系統，并以此為基礎利用保辛攝動的思想，提出一種Hamilton系統的保辛守恒積分算法，實現即保辛又保能量的算法.

圖7為能量的相對誤差，其中計算相對誤差所參考的真實能量H=34.868 J；圖8為約束最大相對誤差.由圖7和8可知，在積分1 000 s的時間內，本文算法計算結果的能量相對誤差和約束相對誤差都相當小，其中能量的相對誤差為10-13數量級，而約束誤差則為10-15數量級，說明本文算法完全是既保約束又保能量的，沒有約束違約問題，計算結果驗證理論分析結果.由于本算例雖然是在空間運動的雙擺，但是給定的初始位移和初始速度設定其實際運動只在一個平面內.圖9和10分別給出質點1的運動軌跡和質點2相對于質點1的運動軌跡.圖9和10中分別給出的3幅圖是從3個不同角度繪制的運動軌跡，其中，圖9（a）為xy平面內的運動軌跡，可見在計算1 000 s的時間區段內，質點1仍然精確地維持在一個平面內運動；同樣圖10（a）亦可驗證這一點.這進一步說明本文方法既保約束又保能量的優點，因為不守恒算法計算時，常常是在最初幾秒可以維持在同一平面；隨著計算繼續，誤差累計，運動軌跡往往躍出平面外.

3結束語

綜合運用文獻[78]提出的祖沖之類方法和文獻[10]提出的保辛守恒思想以及歐拉中點保辛差分格式，建立求解多體動力學方程的既保約束又保能量的算法.該方法在格點處可以嚴格滿足位移約束方程和能量守恒條件，而在非格點處，約束條件和能量方程可以近似滿足，避免約束違約的問題.數值算例表明，本文方法的計算結果滿意.

參考文獻：

[1]胡繼云. 建立多剛體系統動力學方程的坐標變換法及其應用[D]. 重慶： 重慶大學， 2004.

[2]田強. 基于絕對節點坐標方法的柔性多體系統動力學研究與應用[D]. 武漢： 華中科技大學， 2009.

[3]原亮明， 王成國， 劉金朝， 等. 一種求解多體系統微分代數方程的拉格朗日乘子方法[J]. 中國鐵道科學， 2001， 22（2）： 5154.

[4]趙維加， 潘振寬， 王藝兵. 多體系統動力學微分/代數方程約束誤差小擾動自我穩定方法[J]. 應用數學和力學， 2000， 21（1）： 9498.

[5]付士慧， 王琪. 多體系統動力學方程違約修正的數值計算方法[J]. 計算力學學報， 2007， 24（1）： 4449.

[6]戈新生， 趙維加， 陳立群. 基于完全笛卡爾坐標的多體系統微分代數方程符號線性化方法[J]. 工程力學， 2004， 21（4）： 106111.

[7]鐘萬勰， 高強. 約束動力系統的分析結構力學積分[J]. 動力學與控制學報， 2006， 4（3）： 193200.

[8]鐘萬勰， 高強， 彭海軍. 經典力學辛講[M]. 大連： 大連理工大學出版社， 2013.

[9]ZHONG G， MARSDEN J E. LiePoisson HamiltonJacobi theory and LiePoisson integrators[J]. Phys Lett A， 1988， 133（3）： 134139.

[10]高強， 鐘萬勰. Hamilton系統的保辛守恒積分算法[J]. 動力學與控制學報， 2009， 7（3）： 193199.

[11]邢譽峰， 楊蓉. 動力學平衡方程的Euler中點辛差分求解格式[J]. 力學學報， 2007， 39（1）： 100105.

（編輯武曉英）

摘要：針對一類多體動力學問題導出的微分代數方程，提出一種保能量、保約束的算法.該算法基于祖沖之類方法和歐拉中點保辛差分，利用祖沖之類方法保證在時間格點上精確滿足約束方程，避免約束違約問題；并進一步證明該算法在時間格點上可以精確保能量.數值算例進一步驗證該算法的可靠性.

關鍵詞：多體動力學方程； 微分代數方程； 保辛； 祖沖之

中圖分類號： O313.7

文獻標志碼： B

0引言

近年來，多體動力學研究已逐漸成為研究熱點.一方面是由于多體動力學建模得到的往往是強非線性方程，分析困難；另一方面則是因為多體動力學問題具有十分廣泛的運用背景[12]，如航空航天器、車輛和機器人等.目前，多體動力學的建模方法可以分為相對坐標法和絕對坐標法2類.當用絕對坐標法時，得到的往往是微分代數方程組.關于微分代數方程組的求解，目前已有一些進展.

許多學者認為，數值積分過程中約束方程的違約是造成積分困難的重要原因之一.原亮明等[3]把位移約束方程按照泰勒級數展開，與動力學方程組合進行迭代求解；趙維加等[4]用泰勒級數把約束方程展開，根據積分步長提出一種能對約束誤差自動修正的小擾動違約穩定方法；付士慧等[5]對具有完整、定常約束的多體系統，通過修改的帶乘子的拉格朗日正則形式方程，給出一種違約修正方法；戈新生等[6]提出一種基于完全笛卡爾坐標的多體系統微分代數方程符號線性化方法.在諸多研究中，文獻[7]中所提的方法能很好地解決約束方程違約問題.他們利用保辛的時間有限元結合祖沖之類方法的思想：保辛方法可以保證長時間積分的精度，祖沖之類方法可以解決約束方程違約這一問題.但是，當時的文獻中還沒有正式提出祖沖之類方法這一名詞，直到在最近出版的《經典力學辛講》[8]中，祖沖之類方法才被正式提出.祖沖之類方法指出，約束方程不必處處滿足，只要在時間格點處滿足即可.依據這一思想，使得微分代數方程的求解格式簡單，而且無約束違約問題.

對于非線性系統的差分格式要保辛，而國外論文提出，對于不可積系統，保辛和守恒難以同時達成[9]的問題，闡明保辛則能量不能守恒，能量守恒就不能保辛的兩難命題.實際上，文獻[10]通過引入含參變量的近似Hamilton系統，并以此為基礎利用保辛攝動的思想，提出一種Hamilton系統的保辛守恒積分算法，實現即保辛又保能量的算法.

圖7為能量的相對誤差，其中計算相對誤差所參考的真實能量H=34.868 J；圖8為約束最大相對誤差.由圖7和8可知，在積分1 000 s的時間內，本文算法計算結果的能量相對誤差和約束相對誤差都相當小，其中能量的相對誤差為10-13數量級，而約束誤差則為10-15數量級，說明本文算法完全是既保約束又保能量的，沒有約束違約問題，計算結果驗證理論分析結果.由于本算例雖然是在空間運動的雙擺，但是給定的初始位移和初始速度設定其實際運動只在一個平面內.圖9和10分別給出質點1的運動軌跡和質點2相對于質點1的運動軌跡.圖9和10中分別給出的3幅圖是從3個不同角度繪制的運動軌跡，其中，圖9（a）為xy平面內的運動軌跡，可見在計算1 000 s的時間區段內，質點1仍然精確地維持在一個平面內運動；同樣圖10（a）亦可驗證這一點.這進一步說明本文方法既保約束又保能量的優點，因為不守恒算法計算時，常常是在最初幾秒可以維持在同一平面；隨著計算繼續，誤差累計，運動軌跡往往躍出平面外.

3結束語

綜合運用文獻[78]提出的祖沖之類方法和文獻[10]提出的保辛守恒思想以及歐拉中點保辛差分格式，建立求解多體動力學方程的既保約束又保能量的算法.該方法在格點處可以嚴格滿足位移約束方程和能量守恒條件，而在非格點處，約束條件和能量方程可以近似滿足，避免約束違約的問題.數值算例表明，本文方法的計算結果滿意.

參考文獻：

[1]胡繼云. 建立多剛體系統動力學方程的坐標變換法及其應用[D]. 重慶： 重慶大學， 2004.

[2]田強. 基于絕對節點坐標方法的柔性多體系統動力學研究與應用[D]. 武漢： 華中科技大學， 2009.

[3]原亮明， 王成國， 劉金朝， 等. 一種求解多體系統微分代數方程的拉格朗日乘子方法[J]. 中國鐵道科學， 2001， 22（2）： 5154.

[4]趙維加， 潘振寬， 王藝兵. 多體系統動力學微分/代數方程約束誤差小擾動自我穩定方法[J]. 應用數學和力學， 2000， 21（1）： 9498.

[5]付士慧， 王琪. 多體系統動力學方程違約修正的數值計算方法[J]. 計算力學學報， 2007， 24（1）： 4449.

[6]戈新生， 趙維加， 陳立群. 基于完全笛卡爾坐標的多體系統微分代數方程符號線性化方法[J]. 工程力學， 2004， 21（4）： 106111.

[7]鐘萬勰， 高強. 約束動力系統的分析結構力學積分[J]. 動力學與控制學報， 2006， 4（3）： 193200.

[8]鐘萬勰， 高強， 彭海軍. 經典力學辛講[M]. 大連： 大連理工大學出版社， 2013.

[9]ZHONG G， MARSDEN J E. LiePoisson HamiltonJacobi theory and LiePoisson integrators[J]. Phys Lett A， 1988， 133（3）： 134139.

[10]高強， 鐘萬勰. Hamilton系統的保辛守恒積分算法[J]. 動力學與控制學報， 2009， 7（3）： 193199.

[11]邢譽峰， 楊蓉. 動力學平衡方程的Euler中點辛差分求解格式[J]. 力學學報， 2007， 39（1）： 100105.

（編輯武曉英）