การสร้างแบบจำลองวงโคจรคือกระบวนการสร้างแบบจำลองทางคณิตศาสตร์เพื่อจำลองการเคลื่อนที่ของวัตถุขนาดใหญ่ขณะโคจรรอบวัตถุขนาดใหญ่อีกวัตถุหนึ่งเนื่องจากแรงโน้มถ่วง แรงอื่นๆ เช่น แรงดึงดูดจากวัตถุที่สาม แรงต้านอากาศแรงดันจากแสงอาทิตย์หรือแรงขับจาก ระบบ ขับเคลื่อนมักจะถูกจำลองเป็นผลกระทบรอง การสร้างแบบจำลองวงโคจรโดยตรงอาจผลักดันขีดจำกัดความแม่นยำของเครื่องจักรเนื่องจากจำเป็นต้องจำลองการรบกวนเล็กน้อยไปจนถึงวงโคจรขนาดใหญ่มาก ด้วยเหตุนี้ วิธี การรบกวนจึงมักถูกนำมาใช้ในการสร้างแบบจำลองวงโคจรเพื่อให้ได้ความแม่นยำที่ดีขึ้น

การศึกษาการเคลื่อนที่ของวงโคจรและการสร้างแบบจำลองทางคณิตศาสตร์ของวงโคจรเริ่มต้นด้วยความพยายามครั้งแรกในการทำนายการเคลื่อนที่ของดาวเคราะห์บนท้องฟ้า แม้ว่าในสมัยโบราณสาเหตุยังคงเป็นปริศนานิวตันในขณะที่เขากำหนดกฎการเคลื่อนที่และแรงโน้มถ่วง ของเขา ได้นำกฎเหล่านั้นไปใช้ในการวิเคราะห์การรบกวนครั้งแรก^{[ 1 ]}โดยตระหนักถึงความยากลำบากที่ซับซ้อนของการคำนวณ^{[ 1 ]} นักคณิตศาสตร์ผู้ยิ่งใหญ่หลายคนนับตั้งแต่นั้นมาได้ให้ความสนใจกับปัญหาต่างๆ ที่เกี่ยวข้อง ตลอดศตวรรษที่ 18 และ 19 มีความต้องการตารางตำแหน่งของดวงจันทร์และดาวเคราะห์ที่แม่นยำเพื่อวัตถุประสงค์ในการนำทางทางทะเล

การเคลื่อนที่ที่ซับซ้อนของวงโคจรสามารถแยกย่อยได้ การเคลื่อนที่สมมุติที่วัตถุเคลื่อนที่ภายใต้แรงโน้มถ่วงของวัตถุอื่นเพียงวัตถุเดียว โดยทั่วไปจะเป็นภาคตัดกรวยและสามารถจำลองได้ง่ายด้วยวิธีการทางเรขาคณิตนี่เรียกว่าปัญหาของวัตถุสองชิ้นหรือวงโคจรแบบเคปเลอร์ ที่ไม่ถูกรบกวน ความแตกต่างระหว่างวงโคจรแบบเคปเลอร์กับการเคลื่อนที่จริงของวัตถุเกิดจากการรบกวนการรบกวนเหล่านี้เกิดจากแรงอื่นนอกเหนือจากแรงโน้มถ่วงระหว่างวัตถุหลักและวัตถุรอง และต้องจำลองเพื่อสร้างการจำลองวงโคจรที่แม่นยำ วิธีการสร้างแบบจำลองวงโคจรส่วนใหญ่จะจำลองปัญหาของวัตถุสองชิ้นก่อน จากนั้นจึงเพิ่มแบบจำลองของแรงรบกวนเหล่านี้และจำลองแบบจำลองเหล่านี้เมื่อเวลาผ่านไป แรงรบกวนอาจรวมถึงแรงดึงดูดจากวัตถุอื่นนอกเหนือจากวัตถุหลัก ลมสุริยะ แรงต้าน แรงสนามแม่เหล็ก และแรงขับเคลื่อน

วิธีแก้ปัญหาเชิงวิเคราะห์ (นิพจน์ทางคณิตศาสตร์เพื่อทำนายตำแหน่งและการเคลื่อนที่ในอนาคต) สำหรับปัญหาสองและสามวัตถุ แบบง่าย มีอยู่แล้ว แต่ยังไม่พบวิธีแก้ปัญหาสำหรับ ปัญหา nวัตถุยกเว้นในกรณีพิเศษบางกรณี แม้แต่ปัญหาสองวัตถุก็ไม่สามารถแก้ได้หากวัตถุหนึ่งมีรูปร่างไม่สม่ำเสมอ^{[ 2 ]}

เนื่องจากความยากลำบากในการหาคำตอบเชิงวิเคราะห์สำหรับปัญหาที่สนใจส่วนใหญ่การสร้างแบบจำลองและการจำลองด้วย คอมพิวเตอร์ จึงมักถูกนำมาใช้ในการวิเคราะห์การเคลื่อนที่ในวงโคจร มีซอฟต์แวร์หลากหลายประเภทที่สามารถใช้จำลองวงโคจรและวิถีการโคจรของยานอวกาศได้

ในรูปแบบที่ง่ายที่สุด แบบจำลองวงโคจรสามารถสร้างขึ้นได้โดยการสมมติว่ามีเพียงสองวัตถุที่เกี่ยวข้อง ทั้งสองมีพฤติกรรมเหมือนมวลจุดทรงกลม และไม่มีแรงอื่นใดกระทำต่อวัตถุเหล่านั้น ในกรณีนี้ แบบจำลองจะถูกลดทอนให้เหลือเพียงวงโคจรของเคปเลอร์

วงโคจรแบบเคปเลอร์เป็นไปตามภาคตัดกรวยแบบจำลองทางคณิตศาสตร์ของวงโคจรซึ่งแสดงระยะห่างระหว่างวัตถุศูนย์กลางกับวัตถุที่โคจร สามารถแสดงได้ดังนี้:

${\displaystyle r(\nu )={\frac {a(1-e^{2})}{1+e\cos(\nu )}}}$

ที่ไหน:

${\displaystyle r}$คือระยะทาง

${\displaystyle a}$คือแกนกึ่งเอกซึ่งเป็นตัวกำหนดขนาดของวงโคจร

${\displaystyle e}$คือค่าความเยื้องศูนย์กลางซึ่งเป็นตัวกำหนดรูปร่างของวงโคจร

${\displaystyle \nu }$ค่าความผิดปกติที่แท้จริงคือ มุมระหว่างตำแหน่งปัจจุบันของวัตถุที่โคจรอยู่กับตำแหน่งในวงโคจรที่อยู่ใกล้กับวัตถุศูนย์กลางมากที่สุด (เรียกว่า จุดใกล้ที่สุดของวงโคจร )

อีกวิธีหนึ่ง สมการสามารถแสดงได้ดังนี้:

${\displaystyle r(\nu )={\frac {p}{1+e\cos(\nu )}}}$

โดยที่เรียกว่ากึ่งแกนลาตัสเรกตัมของเส้นโค้ง สมการรูปแบบนี้มีประโยชน์อย่างยิ่งเมื่อต้องจัดการกับวิถีโค้งแบบพาราโบลา ซึ่งแกนกึ่งเอกมีค่าเป็นอนันต์ ${\displaystyle p}$

อีกแนวทางหนึ่งใช้กฎแรงโน้มถ่วงสากลของไอแซค นิวตันดังที่นิยามไว้ด้านล่าง:

${\displaystyle F=G{\frac {m_{1}m_{2}}{r^{2}}}}$

ที่ไหน:

${\displaystyle F}$คือขนาดของแรงโน้มถ่วงระหว่างมวลจุดสองจุด

${\displaystyle G}$คือค่าคงที่ความโน้มถ่วง

${\displaystyle m_{1}}$คือมวลของจุดมวลแรก

${\displaystyle m_{2}}$คือมวลของจุดมวลที่สอง

${\displaystyle r}$คือระยะห่างระหว่างจุดมวลสองจุด

หากสมมติเพิ่มเติมว่ามวลของวัตถุหลักมีค่ามากกว่ามวลของวัตถุรองมาก และแทนค่าด้วยกฎการเคลื่อนที่ข้อที่สอง ของนิวตัน จะได้สมการเชิงอนุพันธ์ดังต่อไปนี้

${\displaystyle {\ddot {\mathbf {r} }}={\frac {Gm_{1}}{r^{2}}}\mathbf {\hat {r}} }$

การแก้สม การเชิงอนุพันธ์นี้จะส่งผลให้ได้การเคลื่อนที่แบบเคปเลอร์สำหรับวงโคจร ในทางปฏิบัติ วงโคจรแบบเคปเลอร์มักมีประโยชน์เฉพาะสำหรับการประมาณค่าอันดับแรก กรณีพิเศษ หรือเป็นแบบจำลองพื้นฐานสำหรับวงโคจรที่ถูกรบกวนเท่านั้น

แบบจำลองวงโคจรโดยทั่วไปจะแพร่กระจายในเวลาและอวกาศโดยใช้ วิธี การรบกวน พิเศษ ซึ่งดำเนินการโดยการสร้างแบบจำลองวงโคจรเป็นวงโคจรแบบเคปเลอร์ก่อน จากนั้นจึงเพิ่มการรบกวนเข้าไปในแบบจำลองเพื่ออธิบายการรบกวนต่างๆ ที่ส่งผลต่อวงโคจร^{[ 1 ]} การรบกวนพิเศษสามารถนำไปใช้กับปัญหาใดๆ ในกลศาสตร์ท้องฟ้าได้ เนื่องจากไม่ได้จำกัดเฉพาะกรณีที่แรงรบกวนมีขนาดเล็ก^{[ 2 ]} วิธีการรบกวนพิเศษเป็นพื้นฐานของ ปฏิทินดาวเคราะห์ที่สร้างโดยเครื่องจักรที่มีความแม่นยำที่สุด^{[ 1 ]} ดูตัวอย่างเช่นปฏิทินการพัฒนาของห้องปฏิบัติการเจ็ทโพรพัลชัน

วิธีของ Cowell เป็นวิธีการรบกวนแบบพิเศษ^{[ 3 ]} ในทางคณิตศาสตร์ สำหรับวัตถุที่มีปฏิสัมพันธ์กันแรงนิวตัน บนวัตถุ จากวัตถุอื่นจะถูกรวมเข้าด้วยกันดังนี้ ${\displaystyle n}$${\displaystyle i}$${\displaystyle j}$

${\displaystyle \mathbf {\ddot {r}} _{i}=\sum _{\underset {j\neq i}{j=1}}^{n}{Gm_{j}(\mathbf {r} _{j}-\mathbf {r} _{i}) \over r_{ij}^{3}}}$

ที่ไหน

${\displaystyle \mathbf {\ddot {r}} _{i}}$คือ เวกเตอร์ ความเร่งของวัตถุ${\displaystyle i}$

${\displaystyle G}$คือค่าคงที่ความโน้มถ่วง

${\displaystyle m_{j}}$คือมวลของร่างกาย${\displaystyle j}$

${\displaystyle \mathbf {r} _{i}}$และเป็นเวกเตอร์ตำแหน่งของวัตถุและ${\displaystyle \mathbf {r} _{j}}$${\displaystyle i}$${\displaystyle j}$

${\displaystyle r_{ij}}$คือระยะห่างระหว่างวัตถุกับวัตถุ${\displaystyle i}$${\displaystyle j}$

โดยเวกเตอร์ ทั้งหมด อ้างอิงถึงจุดศูนย์กลางมวลของระบบ สมการนี้จะถูกแยกออกเป็นส่วนประกอบใน, , และส่วนประกอบเหล่านี้จะถูกรวมเข้าด้วยกันในเชิงตัวเลขเพื่อสร้างเวกเตอร์ความเร็วและตำแหน่งใหม่เมื่อการจำลองดำเนินไปข้างหน้าตามเวลา ข้อดีของวิธีการของ Cowell คือความง่ายในการใช้งานและการเขียนโปรแกรม ข้อเสียคือเมื่อการรบกวนมีขนาดใหญ่ (เช่นเมื่อวัตถุเข้าใกล้วัตถุอื่น) ข้อผิดพลาดของวิธีการก็จะใหญ่ขึ้นด้วย^{[ 4 ​​]} ข้อเสียอีกประการหนึ่งคือในระบบที่มีวัตถุกลางที่เด่นชัด เช่นดวงอาทิตย์จำเป็นต้องใช้ตัวเลขสำคัญ จำนวนมาก ในการคำนวณเนื่องจากความแตกต่างอย่างมากของแรงของวัตถุกลางและวัตถุที่รบกวน^{[ 5 ]}${\displaystyle x}$${\displaystyle y}$${\displaystyle z}$

วิธีของ Encke เริ่มต้นด้วยวงโคจรแบบออสคิวเลตเป็นข้อมูลอ้างอิง และทำการอินทิเกรตเชิงตัวเลขเพื่อหาค่าความแปรผันจากข้อมูลอ้างอิงตามฟังก์ชันของเวลา^{[ 6 ]} ข้อดีของวิธีนี้คือโดยทั่วไปแล้วการรบกวนจะมีขนาดเล็ก ดังนั้นการอินทิเกรตจึงสามารถดำเนินการได้ในขั้นตอนที่ใหญ่ขึ้น (ส่งผลให้มีข้อผิดพลาดน้อยลง) และวิธีนี้ได้รับผลกระทบจากการรบกวนที่รุนแรงน้อยกว่าวิธีของ Cowell มาก ข้อเสียคือความซับซ้อน ไม่สามารถใช้งานได้อย่างไม่มีกำหนดโดยไม่ต้องอัปเดตวงโคจรแบบออสคิวเลตเป็นครั้งคราวและดำเนินการต่อจากนั้น ซึ่งเป็นกระบวนการที่เรียกว่าการแก้ไข^{[ 4 ] [ 7 ]}

ให้เป็นเวกเตอร์รัศมีของวงโคจรสัมผัสเป็นเวกเตอร์รัศมีของวงโคจรที่ถูกรบกวน และเป็นค่าเบี่ยงเบนจากวงโคจรสัมผัส ${\displaystyle {\boldsymbol {\rho }}}$${\displaystyle \mathbf {r} }$${\displaystyle \delta \mathbf {r} }$

${\displaystyle \delta \mathbf {r} =\mathbf {r} -{\boldsymbol {\rho }},}$และสมการการเคลื่อนที่ของก็คือ ง่ายๆ เลย${\displaystyle \delta \mathbf {r} }$

1

${\displaystyle {\ddot {\delta \mathbf {r} }}=\mathbf {\ddot {r}} -{\boldsymbol {\ddot {\rho }}}.}$

2

${\displaystyle \mathbf {\ddot {r}} }$และเป็นเพียงสมการการเคลื่อนที่ของและ, ${\displaystyle {\boldsymbol {\ddot {\rho }}}}$${\displaystyle \mathbf {r} }$${\displaystyle {\boldsymbol {\rho }}}$

${\displaystyle \mathbf {\ddot {r}} =\mathbf {a} _{\text{per}}-{\mu \over r^{3}}\mathbf {r} }$สำหรับวงโคจรที่ถูกรบกวนและ

3

${\displaystyle {\boldsymbol {\ddot {\rho }}}=-{\mu \over \rho ^{3}}{\boldsymbol {\rho }}}$สำหรับวงโคจรที่ไม่ถูกรบกวน

4

โดยที่เป็นพารามิเตอร์ความโน้มถ่วงโดยที่และ คือมวลของวัตถุศูนย์กลางและวัตถุที่ถูกรบกวนเป็นความเร่ง ที่ก่อให้ เกิด การรบกวน และและคือขนาดของและ${\displaystyle \mu =G(M+m)}$${\displaystyle M}$${\displaystyle m}$${\displaystyle \mathbf {a} _{\text{per}}}$${\displaystyle r}$${\displaystyle \rho }$${\displaystyle \mathbf {r} }$${\displaystyle {\boldsymbol {\rho }}}$

แทนค่าจากสมการ ( 3 ) และ ( 4 ) ลงในสมการ ( 2 )

${\displaystyle {\ddot {\delta \mathbf {r} }}=\mathbf {a} _{\text{per}}+\mu \left({{\boldsymbol {\rho }} \over \rho ^{3}}-{\mathbf {r} \over r^{3}}\right),}$

5

ซึ่งในทางทฤษฎีสามารถรวมเข้าด้วยกันสองครั้งเพื่อหาค่าเนื่องจากวงโคจรสัมผัสสามารถคำนวณได้ง่ายด้วยวิธีการสองวัตถุและได้รับการพิจารณาและสามารถแก้ไขได้แล้ว ในทางปฏิบัติ ปริมาณในวงเล็บคือผลต่างของเวกเตอร์สองตัวที่เกือบเท่ากัน และจำเป็นต้องมีการจัดการเพิ่มเติมเพื่อหลีกเลี่ยงความจำเป็นในการใช้ตัวเลขสำคัญเพิ่มเติม^{[ 8 ] [ 9 ]}${\displaystyle \delta \mathbf {r} }$${\displaystyle {\boldsymbol {\rho }}}$${\displaystyle \delta \mathbf {r} }$${\displaystyle \mathbf {r} }$${\displaystyle {{\boldsymbol {\rho }} \over \rho ^{3}}-{\mathbf {r} \over r^{3}}}$

ในปี พ.ศ. 2534 Victor R. Bond และ Michael F. Fraietta ได้สร้างวิธีการที่มีประสิทธิภาพและแม่นยำสูงสำหรับการแก้ปัญหาการรบกวนของระบบสองวัตถุ^{[ 10 ]}วิธีนี้ใช้สมการเชิงอนุพันธ์เชิงเส้นและแบบปกติของการเคลื่อนที่ที่ได้มาจาก Hans Sperling และทฤษฎีการรบกวนตามสมการเหล่านี้ที่พัฒนาโดย CA Burdet ในปี พ.ศ. 2407 ในปี พ.ศ. 2516 Bond และ Hanssen ได้ปรับปรุงชุดสมการเชิงอนุพันธ์ของ Burdet โดยใช้พลังงานรวมของระบบที่ถูกรบกวนเป็นพารามิเตอร์แทนพลังงานของระบบสองวัตถุ และลดจำนวนองค์ประกอบเหลือ 13 ในปี พ.ศ. 2532 Bond และ Gottlieb ได้ฝังปริพันธ์ Jacobian ซึ่งเป็นค่าคงที่เมื่อฟังก์ชันศักย์ขึ้นอยู่กับเวลาและตำแหน่งอย่างชัดเจนในสมการของนิวตัน ค่าคงที่ Jacobian ถูกใช้เป็นองค์ประกอบเพื่อแทนที่พลังงานรวมในการกำหนดสมการเชิงอนุพันธ์ของการเคลื่อนที่ใหม่ ในกระบวนการนี้ องค์ประกอบอื่นซึ่งเป็นสัดส่วนกับส่วนประกอบของโมเมนตัมเชิงมุมจะถูกนำมาใช้ ซึ่งทำให้จำนวนองค์ประกอบทั้งหมดกลับมาเป็น 14 ในปี พ.ศ. 2534 บอนด์และเฟรียตตาได้ทำการแก้ไขเพิ่มเติมโดยแทนที่เวกเตอร์ลาปลาสด้วยเวกเตอร์อินทิกรัลอีกตัวหนึ่ง รวมถึงอินทิกรัลสเกลาร์อีกตัวหนึ่งซึ่งกำจัดเทอมฆราวาสขนาดเล็กที่ปรากฏในสมการเชิงอนุพันธ์สำหรับองค์ประกอบบางส่วน^{[ 11 ]}

วิธีการ Sperling–Burdet ดำเนินการตามกระบวนการ 5 ขั้นตอนดังต่อไปนี้: ^{[ 11 ]}

ขั้นตอนที่ 1: การเริ่มต้นใช้งาน

เมื่อกำหนดตำแหน่งเริ่มต้น, ความเร็วเริ่มต้น, และเวลาเริ่มต้น, ตัวแปรต่อไปนี้จะถูกกำหนดค่าเริ่มต้น:${\displaystyle \mathbf {r} _{0}}$${\displaystyle \mathbf {v} _{0}}$${\displaystyle t_{0}}$

${\displaystyle s=0}$

${\displaystyle r_{0}=(\mathbf {r} _{0}\cdot \mathbf {r} _{0})^{1/2}}$

${\displaystyle a=r_{0}}$

${\displaystyle b=\mathbf {r} _{0}\cdot \mathbf {v} _{0}}$

${\displaystyle \tau =t_{0}}$

${\displaystyle {\boldsymbol {\alpha }}=\mathbf {r} _{0}}$

${\displaystyle {\boldsymbol {\beta }}=a\mathbf {v} _{0}}$

การเปลี่ยนแปลงที่เกิดจากมวลรบกวน ซึ่งกำหนดโดยและจะถูกประเมิน${\displaystyle V_{0}}$${\displaystyle {\Bigg [}{\partial {V} \over {\partial {\mathbf {r} }}}{\Bigg ]}_{0}}$

การเปลี่ยนแปลงที่เกิดจากความเร่งอื่นๆ ซึ่งกำหนดโดย จะถูกประเมิน${\displaystyle \mathbf {P} _{0}}$

${\displaystyle \alpha _{J}={\frac {2\mu }{r_{0}}}-\mathbf {v} _{0}\cdot \mathbf {v} _{0}-2V_{0}}$

${\displaystyle \gamma =\mu -\alpha _{J}a}$

${\displaystyle {\boldsymbol {\delta }}=-(\mathbf {v} _{0}\cdot \mathbf {v} _{0})\mathbf {r} _{0}+(\mathbf {r} _{0}\cdot \mathbf {v} _{0})\mathbf {v} _{0}+{\frac {\mu }{r_{0}}}\mathbf {r} _{0}-\alpha _{J}\mathbf {r} _{0}}$

${\displaystyle \sigma =0}$

ขั้นตอนที่ 2: แปลงองค์ประกอบเป็นพิกัด

${\displaystyle \mathbf {r} ={\boldsymbol {\alpha }}+{\boldsymbol {\beta }}sc_{1}+{\boldsymbol {\delta }}s^{2}c_{2}}$

${\displaystyle \mathbf {r'} ={\boldsymbol {\beta }}c_{0}+{\boldsymbol {\delta }}sc_{1}}$

${\displaystyle \mathbf {x} _{3}=\alpha _{J}({\boldsymbol {\alpha }}-\mathbf {r} )+{\boldsymbol {\delta }}}$

${\displaystyle \gamma =\mu -\alpha _{J}a}$

${\displaystyle r=a+bsc_{1}+\gamma s^{2}c_{2}}$

${\displaystyle \mathbf {v} =\mathbf {r'} /r}$

${\displaystyle r'=bc_{0}+\gamma sc_{1}}$

${\displaystyle t=\tau +as+bs^{2}c_{2}+\gamma s^{3}c_{3}}$

ฟังก์ชัน Stumpffอยู่ที่ไหน${\displaystyle c_{0},c_{1},c_{2},c_{3}}$

ขั้นตอนที่ 3: ประเมินสมการเชิงอนุพันธ์สำหรับองค์ประกอบต่างๆ

${\displaystyle \mathbf {F} =\mathbf {P} -{\partial {V} \over \partial {\mathbf {r} }}}$

${\displaystyle \mathbf {Q} =r^{2}\mathbf {F} +2\mathbf {r} (-V+\sigma )}$

${\displaystyle \alpha '_{J}=2(-\mathbf {r'} +r{\boldsymbol {\omega }}\times \mathbf {r} )\cdot \mathbf {P} }$

${\displaystyle \mu {\boldsymbol {\epsilon }}'=2(\mathbf {r'} \cdot \mathbf {F} )\mathbf {r} -(\mathbf {r} \cdot \mathbf {F} )\mathbf {r'} -(\mathbf {r} \cdot \mathbf {r'} )\mathbf {F} }$

${\displaystyle {\boldsymbol {\alpha }}'=-\mathbf {Q} sc_{1}-\mu {\boldsymbol {\epsilon }}'s^{2}c_{2}-\alpha '_{J}{\big [}{\boldsymbol {\alpha }}s^{2}c_{2}+2{\boldsymbol {\beta }}s^{3}{\bar {c}}_{3}+{\frac {1}{2}}{\boldsymbol {\delta }}s^{4}c_{2}^{2}{\big ]}}$

${\displaystyle {\boldsymbol {\beta }}'=\mathbf {Q} c_{0}+\mu {\boldsymbol {\epsilon }}'sc_{1}+\alpha '_{J}{\big [}{\boldsymbol {\alpha }}sc_{1}+{\boldsymbol {\beta }}s^{2}{\bar {c}}_{2}-{\boldsymbol {\delta }}s^{3}(2{\bar {c}}_{3}-c_{1}c_{2}){\big ]}}$

${\displaystyle {\boldsymbol {\delta }}'=\mathbf {Q} \alpha _{J}sc_{1}-\mu {\boldsymbol {\epsilon }}'c_{0}+\alpha '_{J}{\big [}-{\boldsymbol {\alpha }}c_{0}+2\alpha _{J}{\boldsymbol {\beta }}s^{3}{\bar {c}}_{3}+{\frac {1}{2}}{\boldsymbol {\delta }}\alpha _{J}s^{4}c_{2}^{2}{\big ]}}$

${\displaystyle \sigma '=r{\boldsymbol {\omega }}\cdot \mathbf {r} \times \mathbf {F} }$

${\displaystyle a'=-{\frac {1}{r}}\mathbf {r} \cdot \mathbf {Q} sc_{1}-\alpha _{J}'{\big [}as^{2}c_{2}+2bs^{3}{\bar {c}}_{3}+{\frac {1}{2}}\gamma s^{4}c_{2}^{2}{\big ]}}$

${\displaystyle b'={\frac {1}{r}}\mathbf {r} \cdot \mathbf {Q} c_{0}+\alpha _{J}'{\big [}asc_{1}+bs^{2}{\bar {c}}_{2}-\gamma s^{3}(2{\bar {c}}_{3}-c_{1}c_{2}){\big ]}}$

${\displaystyle \gamma '=-{\frac {1}{r}}\mathbf {r} \cdot \mathbf {Q} \alpha _{J}sc_{1}+\alpha _{J}'{\big [}-ac_{0}+2b\alpha _{J}s^{3}{\bar {c}}_{3}+{\frac {1}{2}}\gamma \alpha _{J}s^{4}c_{2}^{2}{\big ]}}$

${\displaystyle \tau '={\frac {1}{r}}\mathbf {r} \cdot \mathbf {Q} s^{2}c_{2}+\alpha _{J}'{\big [}as^{3}c_{3}+{\frac {1}{2}}bs^{4}c_{2}^{2}-2\gamma s^{5}(c_{5}-4{\bar {c}}_{5}){\big ]}}$

ขั้นตอนที่ 4: การบูรณาการ

ในที่นี้ สมการเชิงอนุพันธ์จะถูกอินทิเกรตในช่วงเวลาหนึ่งเพื่อให้ได้ค่าองค์ประกอบที่${\displaystyle \Delta s}$${\displaystyle s+\Delta s}$

ขั้นตอนที่ 5: เลื่อนไปข้างหน้า

ตั้งค่าและกลับไปที่ขั้นตอนที่ 2 จนกว่าจะตรงตามเงื่อนไขการหยุดการจำลอง${\displaystyle s=s+\Delta s}$

แรงรบกวนต่างๆ ทำให้วงโคจรเบี่ยงเบนไปจากวงโคจรแบบเคปเลอร์ที่สมบูรณ์แบบ มีการสร้างและเรียกใช้แบบจำลองสำหรับแรงแต่ละชนิดในระหว่างการจำลองวงโคจร เพื่อให้สามารถกำหนดผลกระทบของแรงเหล่านั้นต่อวงโคจรได้

โลกไม่ได้เป็นทรงกลมที่สมบูรณ์แบบ และมวลก็ไม่ได้กระจายอย่างสม่ำเสมอภายในโลก ส่งผลให้แบบจำลองแรงโน้มถ่วงแบบจุดมวลไม่แม่นยำสำหรับวงโคจรรอบโลก โดยเฉพาะวงโคจรต่ำของโลกเพื่ออธิบายความแปรผันของศักยภาพแรงโน้มถ่วงรอบพื้นผิวโลก สนามแรงโน้มถ่วงของโลกจึงถูกจำลองด้วยฮาร์มอนิกทรงกลม^{[ 12 ]}ซึ่งแสดงผ่านสมการ:

${\displaystyle {\mathbf {f} }=-{\frac {\mu }{R^{2}}}\mathbf {\hat {r}} +\sum _{n=2}^{\infty }\sum _{m=0}^{n}{\mathbf {f} }_{n,m}}$

ที่ไหน

${\displaystyle {\mu }}$คือพารามิเตอร์ความโน้มถ่วง ซึ่งนิยามว่าเป็นผลคูณของ G ซึ่งเป็นค่าคงที่ความโน้มถ่วงสากลและมวลของวัตถุหลัก

${\displaystyle \mathbf {\hat {r}} }$คือเวกเตอร์หน่วยที่กำหนดระยะห่างระหว่างวัตถุหลักและวัตถุรอง โดยที่คือขนาดของระยะห่าง${\displaystyle {R}}$

${\displaystyle {\mathbf {f} }_{n,m}}$แสดงถึงการมีส่วนร่วมของฮาร์มอนิกทรงกลมระดับnและลำดับmซึ่งกำหนดไว้ดังนี้: ^{[ 12 ]}${\displaystyle {\mathbf {f} }}$

${\displaystyle {\begin{aligned}\mathbf {f} _{n,m}&={\frac {\mu R_{O}^{2}}{R^{n+m+1}}}\left({\frac {C_{n,m}{\mathcal {C}}_{m}+S_{n,m}{\mathcal {S}}_{m}}{R}}(A_{n,m+1}\mathbf {\hat {e}} _{3}-\left(s_{\lambda }A_{n,m+1}+(n+m+1)A_{n,m}\right)\mathbf {\hat {r}} \right)\\[10pt]&{}\quad {}+mA_{n,m}((C_{n,m}{\mathcal {C}}_{m-1}+S_{n,m}{\mathcal {S}}_{m-1})\mathbf {\hat {e}} _{1}+(S_{n,m}{\mathcal {C}}_{m-1}-C_{n,m}{\mathcal {S}}_{m-1})\mathbf {\hat {e}} _{2}))\end{aligned}}}$

ที่ไหน:

${\displaystyle R_{O}}$คือรัศมีเฉลี่ยบริเวณเส้นศูนย์สูตรของวัตถุหลัก

${\displaystyle R}$คือขนาดของเวกเตอร์ตำแหน่งจากจุดศูนย์กลางของวัตถุหลักไปยังจุดศูนย์กลางของวัตถุรอง

${\displaystyle C_{n,m}}$และเป็นสัมประสิทธิ์ความโน้มถ่วงระดับnและอันดับmซึ่งโดยทั่วไปจะหาได้จากการวัดความโน้มถ่วง${\displaystyle S_{n,m}}$

เวกเตอร์หน่วยกำหนดระบบพิกัดที่ตรึงอยู่กับวัตถุหลัก สำหรับโลกนั้น เวกเตอร์อยู่ในระนาบเส้นศูนย์สูตร ขนานกับเส้นที่ตัดผ่านจุดศูนย์กลางทางเรขาคณิตของโลกและเส้นเมริเดียนกรีนิช เวกเตอร์ชี้ไปในทิศทางของแกนขั้วโลกเหนือ และ${\displaystyle \mathbf {\hat {e}} _{1},\mathbf {\hat {e}} _{2},\mathbf {\hat {e}} _{3}}$${\displaystyle \mathbf {\hat {e}} _{1}}$${\displaystyle \mathbf {\hat {e}} _{3}}$${\displaystyle \mathbf {\hat {e}} _{2}=\mathbf {\hat {e}} _{3}\times \mathbf {\hat {e}} _{1}}$

${\displaystyle A_{n,m}}$เรียกว่าพหุนามเลอจองเดอร์อนุพันธ์ที่มีดีกรีnและอันดับmสามารถแก้ได้โดยใช้ความสัมพันธ์เวียนเกิดดังนี้:${\displaystyle A_{n,m}(u)={\frac {1}{n-m}}((2n-1)uA_{n-1,m}(u)-(n+m-1)A_{n-2,m}(u))}$

${\displaystyle s_{\lambda }}$คือค่าไซน์ของละติจูดทางภูมิศาสตร์ของวัตถุรอง ซึ่งก็คือ.${\displaystyle \mathbf {\hat {r}} \cdot \mathbf {\hat {e}} _{3}}$

${\displaystyle {\mathcal {C}}_{m},{\mathcal {S}}_{m}}$ถูกกำหนดด้วยความสัมพันธ์เวียนเกิดและเงื่อนไขเริ่มต้นดังต่อไปนี้:${\displaystyle {\mathcal {C}}_{m}={\mathcal {C}}_{1}{\mathcal {C}}_{m-1}-{\mathcal {S}}_{1}{\mathcal {S}}_{m-1},{\mathcal {S}}_{m}={\mathcal {S}}_{1}{\mathcal {C}}_{m-1}+{\mathcal {C}}_{1}{\mathcal {S}}_{m-1},{\mathcal {S}}_{0}=0,{\mathcal {S}}_{1}=\mathbf {R} \cdot \mathbf {\hat {e}} _{2},{\mathcal {C}}_{0}=1,{\mathcal {C}}_{1}=\mathbf {R} \cdot \mathbf {\hat {e}} _{1}}$

เมื่อจำลองการรบกวนของวงโคจรรอบวัตถุหลักจะต้องรวมเฉพาะผลรวมของเทอมต่างๆ ในการรบกวนเท่านั้น เนื่องจากแบบจำลองแรงโน้มถ่วงของมวลจุดได้ถูกนำมาพิจารณาในเทอม นั้นแล้ว${\displaystyle {\mathbf {f} }_{n,m}}$${\displaystyle -{\frac {\mu }{R^{2}}}\mathbf {\hat {r}} }$

แรงโน้มถ่วงจากวัตถุที่สามสามารถทำให้เกิดการรบกวนต่อวงโคจรได้ ตัวอย่างเช่นดวงอาทิตย์และดวงจันทร์ทำให้เกิดการรบกวนต่อวงโคจรรอบโลก^{[ 13 ]}แรงเหล่านี้ถูกจำลองในลักษณะเดียวกับที่แรงโน้มถ่วงถูกจำลองสำหรับวัตถุหลักโดยใช้การจำลองแรงโน้มถ่วง N-body โดยตรงโดยทั่วไปแล้ว จะใช้แบบจำลองแรงโน้มถ่วงแบบจุดมวลทรงกลมเท่านั้นในการจำลองผลกระทบจากวัตถุที่สามเหล่านี้^{[ 14 ]} กรณีพิเศษบางกรณีของการรบกวนจากวัตถุที่สามมีคำตอบเชิงวิเคราะห์โดยประมาณ ตัวอย่างเช่น การรบกวนสำหรับไรต์แอสเซนชันของโหนดขึ้นและอาร์กิวเมนต์ของเพริจีสำหรับวงโคจรโลกแบบวงกลมคือ: ^{[ 13 ]}

${\displaystyle {\dot {\Omega }}_{\mathrm {MOON} }=-0.00338(\cos(i))/n}$

${\displaystyle {\dot {\omega }}_{\mathrm {MOON} }=-0.00169(4-5\sin ^{2}(i))/n}$

ที่ไหน:

${\displaystyle {\dot {\Omega }}}$คือการเปลี่ยนแปลงของไรต์แอสเซนชันของจุดแอสเซนชันในหน่วยองศาต่อวัน

${\displaystyle {\dot {\omega }}}$คือการเปลี่ยนแปลงของมุมของจุดใกล้โลกที่สุดในหน่วยองศาต่อวัน

${\displaystyle i}$คือ ค่าความเอียง ของวงโคจร

${\displaystyle n}$คือจำนวนรอบการโคจรต่อวัน

แรงดันรังสีจากดวงอาทิตย์ทำให้เกิดการรบกวนวงโคจร ขนาดของความเร่งที่ส่งผลต่อยานอวกาศในวงโคจรของโลกนั้นจำลองโดยใช้สมการด้านล่าง: ^{[ 13 ]}

${\displaystyle a_{R}\approx -4.5\times 10^{-6}(1+r)A/m}$

ที่ไหน:

${\displaystyle a_{R}}$คือขนาดของความเร่งในหน่วยเมตรต่อวินาที²

${\displaystyle A}$คือพื้นที่หน้าตัดที่สัมผัสกับแสงแดดในหน่วยตารางเมตร

${\displaystyle m}$มวลของยานอวกาศมีหน่วยเป็นกิโลกรัมเท่าใด

${\displaystyle r}$คือค่าตัวประกอบการสะท้อนซึ่งขึ้นอยู่กับคุณสมบัติของวัสดุสำหรับการดูดกลืนสำหรับการสะท้อนแบบกระจกเงาและสำหรับ การ สะท้อนแบบกระจาย${\displaystyle r=0}$${\displaystyle r=1}$${\displaystyle r\approx 0.4}$

สำหรับวงโคจรรอบโลก แรงดันรังสีจากดวงอาทิตย์จะกลายเป็นแรงที่แข็งแกร่งกว่าแรงต้านเมื่ออยู่เหนือระดับความสูง 800 กม. (500 ไมล์) ^{[ 13 ]}

มีระบบขับเคลื่อนยานอวกาศหลายประเภท เครื่องยนต์จรวดเป็นหนึ่งในระบบที่ใช้กันอย่างแพร่หลาย แรงของเครื่องยนต์จรวดถูกจำลองโดยสมการ: ^{[ 15 ]}

${\displaystyle F_{n}={\dot {m}}\;v_{\text{e}}={\dot {m}}\;v_{\text{e-act}}+A_{\text{e}}(p_{\text{e}}-p_{\text{amb}})}$

ที่ไหน:
${\displaystyle {\dot {m}}}$	= อัตราการไหลของมวลก๊าซไอเสีย
${\displaystyle v_{\text{e}}}$	= ความเร็วไอเสียที่มีประสิทธิภาพ
${\displaystyle v_{\text{e-act}}}$	= ความเร็วไอพ่นจริงที่ระนาบทางออกของหัวฉีด
${\displaystyle A_{\text{e}}}$	= พื้นที่หน้าตัดการไหล ณ ระนาบทางออกของหัวฉีด (หรือระนาบที่ไอพ่นออกจากหัวฉีดหากเป็นการไหลแบบแยกส่วน)
${\displaystyle p_{\text{e}}}$	= ความดันสถิตที่ระนาบทางออกของหัวฉีด
${\displaystyle p_{\text{amb}}}$	= ความดันบรรยากาศ (หรือความดันแวดล้อม)

อีกวิธีหนึ่งที่เป็นไปได้คือใบเรือพลังงานแสงอาทิตย์ใบเรือพลังงานแสงอาทิตย์ใช้แรงดันรังสีในลักษณะที่ทำให้เกิดแรงขับเคลื่อนที่ต้องการ^{[ 16 ]}แบบจำลองการรบกวนเนื่องจากลมสุริยะสามารถใช้เป็นแบบจำลองของแรงขับเคลื่อนจากใบเรือพลังงานแสงอาทิตย์ได้

แรงที่ไม่ใช่แรงโน้มถ่วงหลักที่กระทำต่อดาวเทียมในวงโคจรต่ำของโลกคือแรงต้านของบรรยากาศ^{[ 13 ]}แรงต้านจะกระทำในทิศทางตรงข้ามกับความเร็วและดึงพลังงานออกจากวงโคจร แรงเนื่องจากแรงต้านถูกจำลองโดยสมการต่อไปนี้:

${\displaystyle F_{D}\,=\,{\tfrac {1}{2}}\,\rho \,v^{2}\,C_{d}\,A,}$

ที่ไหน

${\displaystyle \mathbf {F} _{D}}$คือแรงต้าน

${\displaystyle \mathbf {} \rho }$คือความหนาแน่นของของเหลว^{[ a ]}

${\displaystyle \mathbf {} v}$คือความเร็วของวัตถุเมื่อเทียบกับของเหลว

${\displaystyle \mathbf {} C_{d}}$คือค่าสัมประสิทธิ์แรงต้าน ( พารามิเตอร์ไร้มิติ เช่น 2 ถึง 4 สำหรับดาวเทียมส่วนใหญ่^{[ 13 ]} )

${\displaystyle \mathbf {} A}$คือพื้นที่อ้างอิง

โดยทั่วไปแล้ว วงโคจรที่มีระดับความสูงต่ำกว่า 120 กม. (75 ไมล์) จะมีแรงต้านสูงมากจนวงโคจรเสื่อมลงเร็วเกินไป ทำให้ดาวเทียมมีอายุการใช้งานไม่เพียงพอต่อการปฏิบัติภารกิจใดๆ ในทางปฏิบัติ ในทางกลับกัน วงโคจรที่มีระดับความสูงมากกว่า 600 กม. (370 ไมล์) จะมีแรงต้านค่อนข้างน้อย ทำให้วงโคจรเสื่อมลงช้าพอที่จะไม่มีผลกระทบต่อดาวเทียมตลอดอายุการใช้งาน^{[ 13 ]}ความหนาแน่นของอากาศในชั้นเทอร์โมสเฟียร์ซึ่งเป็นที่อยู่ของดาวเทียมที่โคจรรอบโลกในระดับต่ำส่วนใหญ่สามารถเปลี่ยนแปลงได้อย่างมาก การเปลี่ยนแปลงนี้เกิดจากกิจกรรมของดวงอาทิตย์เป็นหลัก ดังนั้นกิจกรรมของดวงอาทิตย์จึงสามารถส่งผลกระทบอย่างมากต่อแรงต้านของยานอวกาศและทำให้การจำลองวงโคจรในระยะยาวมีความซับซ้อนมากขึ้น^{[ 13 ]}

สนามแม่เหล็กสามารถมีบทบาทสำคัญในฐานะแหล่งที่มาของการรบกวนวงโคจร ดังที่เห็นได้ในLong Duration Exposure Facility [ ^{12 ] เช่น}เดียวกับแรงโน้มถ่วง สนามแม่เหล็กของโลกสามารถแสดงได้ผ่านฮาร์มอนิกทรงกลมดังที่แสดงด้านล่าง: ^{[ 12 ]}

${\displaystyle {\mathbf {B} }=\sum _{n=1}^{\infty }\sum _{m=0}^{n}{\mathbf {B} }_{n,m}}$

ที่ไหน

${\displaystyle {\mathbf {B} }}$คือเวกเตอร์สนามแม่เหล็ก ณ จุดหนึ่งเหนือพื้นผิวโลก

${\displaystyle {\mathbf {B} }_{n,m}}$แสดงถึงการมีส่วนร่วมของฮาร์มอนิกทรงกลมระดับnและลำดับmซึ่งกำหนดไว้ดังนี้: ^{[ 12 ]}${\displaystyle {\mathbf {B} }}$

${\displaystyle {\begin{aligned}\mathbf {B} _{n,m}={}&{\frac {K_{n,m}a^{n+2}}{R^{n+m+1}}}\left[{\frac {g_{n,m}{\mathcal {C}}_{m}+h_{n,m}{\mathcal {S}}_{m}}{R}}((s_{\lambda }A_{n,m+1}+(n+m+1)A_{n,m})\mathbf {\hat {r}} )-A_{n,m+1}\mathbf {\hat {e}} _{3}\right]\\[10pt]&{}-mA_{n,m}((g_{n,m}{\mathcal {C}}_{m-1}+h_{n,m}{\mathcal {S}}_{m-1})\mathbf {\hat {e}} _{1}+(h_{n,m}{\mathcal {C}}_{m-1}-g_{n,m}{\mathcal {S}}_{m-1})\mathbf {\hat {e}} _{2}))\end{aligned}}}$

ที่ไหน:

${\displaystyle a}$คือรัศมีเฉลี่ยบริเวณเส้นศูนย์สูตรของวัตถุหลัก

${\displaystyle R}$คือขนาดของเวกเตอร์ตำแหน่งจากจุดศูนย์กลางของวัตถุหลักไปยังจุดศูนย์กลางของวัตถุรอง

${\displaystyle \mathbf {\hat {r}} }$เป็นเวกเตอร์หน่วยที่ชี้ไปในทิศทางของวัตถุรอง โดยมีจุดกำเนิดอยู่ที่ศูนย์กลางของวัตถุหลัก

${\displaystyle g_{n,m}}$และเป็นสัมประสิทธิ์เกาส์ระดับnและอันดับmซึ่งโดยทั่วไปจะหาได้จากการวัดสนามแม่เหล็ก${\displaystyle h_{n,m}}$

เวกเตอร์หน่วยกำหนดระบบพิกัดที่ตรึงอยู่กับวัตถุหลัก สำหรับโลกนั้น เวกเตอร์อยู่ในระนาบเส้นศูนย์สูตร ขนานกับเส้นที่ตัดผ่านจุดศูนย์กลางทางเรขาคณิตของโลกและเส้นเมริเดียนกรีนิช เวกเตอร์ชี้ไปในทิศทางของแกนขั้วโลกเหนือ และ${\displaystyle \mathbf {\hat {e}} _{1},\mathbf {\hat {e}} _{2},\mathbf {\hat {e}} _{3}}$${\displaystyle \mathbf {\hat {e}} _{1}}$${\displaystyle \mathbf {\hat {e}} _{3}}$${\displaystyle \mathbf {\hat {e}} _{2}=\mathbf {\hat {e}} _{3}\times \mathbf {\hat {e}} _{1}}$

${\displaystyle A_{n,m}}$เรียกว่าพหุนามเลอจองเดอร์อนุพันธ์ที่มีดีกรีnและอันดับmสามารถแก้ได้โดยใช้ความสัมพันธ์เวียนเกิดดังนี้:${\displaystyle A_{n,m}(u)={\frac {1}{n-m}}((2n-1)uA_{n-1,m}(u)-(n+m-1)A_{n-2,m}(u))}$

${\displaystyle K_{n,m}}$นิยามว่า: 1 ถ้าm  = 0 สำหรับและและสำหรับและ${\displaystyle {\big [}{\frac {n-m}{n+m}}{\big ]}^{0.5}K_{n-1,m}}$${\displaystyle n\geq (m+1)}$${\displaystyle m=[1\ldots \infty ]}$${\displaystyle [(n+m)(n-m+1)]^{-0.5}K_{n,m-1}}$${\displaystyle n\geq m}$${\displaystyle m=[2\ldots \infty ]}$

${\displaystyle s_{\lambda }}$คือค่าไซน์ของละติจูดทางภูมิศาสตร์ของวัตถุรอง ซึ่งก็คือ.${\displaystyle \mathbf {\hat {r}} \cdot \mathbf {\hat {e}} _{3}}$

${\displaystyle {\mathcal {C}}_{m},{\mathcal {S}}_{m}}$ถูกกำหนดด้วยความสัมพันธ์เวียนเกิดและเงื่อนไขเริ่มต้นดังต่อไปนี้:${\displaystyle {\mathcal {C}}_{m}={\mathcal {C}}_{1}{\mathcal {C}}_{m-1}-{\mathcal {S}}_{1}{\mathcal {S}}_{m-1},{\mathcal {S}}_{m}={\mathcal {S}}_{1}{\mathcal {C}}_{m-1}+{\mathcal {C}}_{1}{\mathcal {S}}_{m-1},{\mathcal {S}}_{0}=0,{\mathcal {S}}_{1}=\mathbf {R} \cdot \mathbf {\hat {e}} _{2},{\mathcal {C}}_{0}=1,{\mathcal {C}}_{1}=\mathbf {R} \cdot \mathbf {\hat {e}} _{1}}$

• ปัญหา n-body
• การสั่นพ้องของวงโคจร
• วงโคจรการสั่น
• การรบกวน (ดาราศาสตร์)
• ขอบเขตอิทธิพล (ดาราศาสตร์พลศาสตร์)
• ปัญหาวัตถุสองชิ้น

• [1]แผนที่แรงโน้มถ่วงของโลก