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

吳鋒　高強　鐘萬勰

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

關(guān)鍵詞：多體動力學(xué)方程； 微分代數(shù)方程； 保辛； 祖沖之

中圖分類號： O313.7

文獻標(biāo)志碼： B

0引言

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

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

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

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

3結(jié)束語

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

參考文獻：

[1]胡繼云. 建立多剛體系統(tǒng)動力學(xué)方程的坐標(biāo)變換法及其應(yīng)用[D]. 重慶： 重慶大學(xué)， 2004.

[2]田強. 基于絕對節(jié)點坐標(biāo)方法的柔性多體系統(tǒng)動力學(xué)研究與應(yīng)用[D]. 武漢： 華中科技大學(xué)， 2009.

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

[4]趙維加， 潘振寬， 王藝兵. 多體系統(tǒng)動力學(xué)微分/代數(shù)方程約束誤差小擾動自我穩(wěn)定方法[J]. 應(yīng)用數(shù)學(xué)和力學(xué)， 2000， 21（1）： 9498.

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

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

[7]鐘萬勰， 高強. 約束動力系統(tǒng)的分析結(jié)構(gòu)力學(xué)積分[J]. 動力學(xué)與控制學(xué)報， 2006， 4（3）： 193200.

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

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

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

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

（編輯武曉英）

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

關(guān)鍵詞：多體動力學(xué)方程； 微分代數(shù)方程； 保辛； 祖沖之

中圖分類號： O313.7

文獻標(biāo)志碼： B

0引言

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

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

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

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

3結(jié)束語

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

參考文獻：

[1]胡繼云. 建立多剛體系統(tǒng)動力學(xué)方程的坐標(biāo)變換法及其應(yīng)用[D]. 重慶： 重慶大學(xué)， 2004.

[2]田強. 基于絕對節(jié)點坐標(biāo)方法的柔性多體系統(tǒng)動力學(xué)研究與應(yīng)用[D]. 武漢： 華中科技大學(xué)， 2009.

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

[4]趙維加， 潘振寬， 王藝兵. 多體系統(tǒng)動力學(xué)微分/代數(shù)方程約束誤差小擾動自我穩(wěn)定方法[J]. 應(yīng)用數(shù)學(xué)和力學(xué)， 2000， 21（1）： 9498.

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

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

[7]鐘萬勰， 高強. 約束動力系統(tǒng)的分析結(jié)構(gòu)力學(xué)積分[J]. 動力學(xué)與控制學(xué)報， 2006， 4（3）： 193200.

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

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

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

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

（編輯武曉英）

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

關(guān)鍵詞：多體動力學(xué)方程； 微分代數(shù)方程； 保辛； 祖沖之

中圖分類號： O313.7

文獻標(biāo)志碼： B

0引言

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

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

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

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

3結(jié)束語

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

參考文獻：

[1]胡繼云. 建立多剛體系統(tǒng)動力學(xué)方程的坐標(biāo)變換法及其應(yīng)用[D]. 重慶： 重慶大學(xué)， 2004.

[2]田強. 基于絕對節(jié)點坐標(biāo)方法的柔性多體系統(tǒng)動力學(xué)研究與應(yīng)用[D]. 武漢： 華中科技大學(xué)， 2009.

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

[4]趙維加， 潘振寬， 王藝兵. 多體系統(tǒng)動力學(xué)微分/代數(shù)方程約束誤差小擾動自我穩(wěn)定方法[J]. 應(yīng)用數(shù)學(xué)和力學(xué)， 2000， 21（1）： 9498.

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

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

[7]鐘萬勰， 高強. 約束動力系統(tǒng)的分析結(jié)構(gòu)力學(xué)積分[J]. 動力學(xué)與控制學(xué)報， 2006， 4（3）： 193200.

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

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

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

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

（編輯武曉英）