صياغة قوة الضغط لتجمعات الغبار من عوامل ملء الحجم المنخفضة إلى العالية باستخدام المحاكاة العددية

Misako Tatsuuma وَ Akimasa Kataoka وَ Satoshi Okuzumi وَ Hidekazu Tanaka

LaTeX

ملخص

تعد قوة الضغط مفتاحًا لفهم البنية الداخلية لتجمعات الغبار في أقراص الكواكب الأولية والأجسام الناتجة عنها، مثل المذنبات والكويكبات في النظام الشمسي. نمذَّجت الأعمال السابقة قوة الضغط لتجمعات الغبار شديدة المسامية عند عوامل ملء حجم أقل من 0.1، ومع ذلك لا يزال الفهم الشامل لقوة الضغط عبر نطاق عوامل ملء الحجم من المنخفضة (\(<0.1\)) إلى العالية (\(>0.1\)) غير مكتمل. في هذه الورقة نجري محاكاة عددية لضغط التجمعات استنادًا إلى نظرية JKR لصياغة قوة الضغط بشكل شامل. نجري سلسلة من المحاكاة العددية باستخدام حدود دورية متحركة لمتابعة سلوك الضغط. من النتائج نجد أن قوة الضغط تزيد شدتها بشكل حاد عند تجاوز عامل ملء الحجم 0.1. نقدم صيغةً شاملة تأخذ في الاعتبار حركة التدحرج للتجمعات عند عوامل الملء المنخفضة والترتيب الأقرب للتجمعات عند عوامل الملء العالية. كما نكتشف أن آليات الضغط السائدة تتغير من التدحرج عند العوامل المنخفضة إلى الانزلاق والالتواء عند العوامل العالية. تتوافق نتائجنا جيدًا مع الدراسات العددية السابقة، وتقترح معادلتنا التحليلية توافقًا مع النتائج التجريبية إذا افترضنا أن طاقة سطح السيليكات ≃210±90 mJ m⁻². يمكن تطبيق نتائجنا الآن على خصائص الأجسام الصلبة الصغيرة مثل المذنبات والكويكبات والحصى.

مقدمة

الخطوة الأولى في تكوّن الكواكب هي تجمع حبيبات الغبار التي يقل حجمها عن الميكرون. يُطلق على هذه البنى اسم "تجمعات الغبار" (Smirnov1990, Meakin1991, Ossenkopf1993, Dominik1997, Wurm1998, Kempf1999, Blum2000, Krause2004, Paszun2006, Paszun2008, Paszun2009, Wada2007, Wada2008, Wada2009, Wada2013, Suyama2008, Suyama2012, Okuzumi2009dustagg, Geretshauser2010, Geretshauser2011). في المرحلة الأولية من نمو الغبار، يصطدم تجمع غبار بآخر ويلتصق به، منتجًا تجمعات كتلية كروية (Mukai1992). تؤدي هذه التجمعات لاحقًا إلى تكوّن كويكبات صغيرة، وهي هياكل يصل حجمها إلى الكيلومترات (Okuzumi2012, Kataoka2013L). هناك سيناريو آخر ينمو فيه الغبار إلى حصى بحجم المليمتر، ثم تتجمع هذه الحصى لتكوين كويكبات صغيرة عبر آليات عدم الاستقرار أو الاصطدام (Johansen2007, Windmark2012b, Davidsson2016, WahlbergJansson2017, Yang2017, Lorek2018, Fulle2019). في هذا السيناريو تختلف البنية الداخلية لتجمعات الحصى عن تلك لتجمعات الغبار التي ندرسها هنا.

الضغط على تجمعات الغبار هو عملية رئيسية أثناء نموها، ويشمل عدة آليات: الاصطدام، ضغط الغاز في القرص، والضغط الناتج عن الجاذبية الذاتية. أظهرت الدراسات العددية أن الضغط الاصطدامي لا يكفي، فتظل الكثافات الداخلية منخفضة حوالي \(\sim10^{-5}\mathrm{\ g\ cm^{-3}}\) (Okuzumi2012, Kataoka2013L). تشير الدراسات التجريبية إلى أن ارتداد التجمعات قد يؤدي إلى الضغط (Krijt2018)، لكن المحاكاة العددية تظهر أن ارتداد التجمعات المسامية نادر (Wada2011). بالنسبة لضغط الغاز والجاذبية الذاتية، تحدد قوة الضغط الكثافة الداخلية للتجمعات (Blum2004, Paszun2008, Guttler2009, Seizinger2012, Kataoka2013, Kataoka2013L, Omura2017)، وتمتد أهميتها أيضًا للأجسام الأكبر مثل الكويكبات والمذنبات (Omura2018, Omura2021).

قدّم (Kataoka2013) نموذجًا تحليليًا لقوة الضغط لتجمعات شديدة المسامية عند عوامل ملء حجمية أقل من 0.1، مستخدمًا عامل التعبئة وعدة معاملات مادية مثل نصف قطر الحبيبات وطاقتها السطحية.

ومع ذلك، لا يزال الفهم الشامل لقوة الضغط من عوامل الملء الحجمية المنخفضة (\(<0.1\)) إلى العالية (\(>0.1\)) مفقودًا. فالقوة عند عوامل الملء العالية ضرورية لتطبيقات المذنبات والكويكبات والحصى، بينما عوامل الملء المنخفضة مهمة لنمو الغبار. تناولت بعض الدراسات قوة الضغط عند عوامل ملء حجمية أعلى من 0.1 (Blum2004, Paszun2008, Guttler2009, Seizinger2012, Omura2017, Omura2018)، لكن الاعتماد على المعاملات المادية لا يزال غير واضح ويظهر فجوة بين النطاقين المنخفض والعالي.

في هذا العمل نجري محاكاة عددية لضغط تجمعات الغبار ونصوغ قوة الضغط التي تغطي كامل نطاق عوامل التعبئة الحجمية. نستخدم نفس الشفرة التي استُخدمت في (Kataoka2013)، لكننا نركز على عوامل الملء العالية (>0.1) لتطبيق النتائج على الأجسام الصغيرة في النظام الشمسي. ندرس أيضًا تبعية القوة على معاملات مادية مثل نصف القطر وطاقة السطح. وأخيرًا نستخلص صيغةً تحليلية مصححة (corrected) يمكن تطبيقها على معاملات مادية أخرى.

ينظم هذا البحث على النحو التالي. في القسم [sec:setting] نقدم إعدادات المحاكاة ونموذج التفاعل الأحادي بناءً على (Dominik1997) و(Wada2007). نشرح شروط الحدود الدورية المتحركة وحساب القوة في [subsec:setting:boundary]، وطريقة قياس عامل التعبئة الحجمي في [subsec:setting:measure]. في القسم [sec:result] نعرض النتائج العددية وتبعياتها على المعاملات. في القسم [sec:discuss] نصحح الصيغة التحليلية ونناقش الفيزياء خلف الضغط وآليات تبديد الطاقة، ثم نقارن نتائجنا مع الدراسات السابقة. نختتم في القسم [sec:conclusion].

إعدادات المحاكاة

في هذا القسم نشرح إعدادات المحاكاة. أولًا، نعرض نموذج التفاعل الأحادي استنادًا إلى (Dominik1997) و(Wada2007) في [subsec:setting:model]، مع توضيح التخميد الاصطناعي. ثانيًا، نصف تصميم المحاكاة الذي يستخدم حدودًا دورية متحركة لضغط التجمعات في [subsec:setting:boundary]، بما في ذلك الشروط الأولية وسرعة الحدود. ثالثًا، نشرح كيفية حساب قوة الضغط وعامل التعبئة الحجمي في [subsec:setting:measure].

نموذج تفاعل الوحدات الأحادية

نحسب تفاعلات الوحدات الأحادية الكروية المتلامسة باستخدام نموذج يستند إلى نظرية (Johnson1971) وتطويرات (Dominik1997, Wada2007). يتضمن النموذج أربعة أنواع من التفاعلات: الانضغاط العادي، الانزلاق، الدوران، والالتواء. المعاملات المادية هي نصف قطر الوحدة الأحادية \(r_0\)، الكثافة \(\rho_0\)، طاقة السطح \(\gamma\)، معامل بواسون \(\nu\)، معامل يونغ \(E\)، والإزاحة الدورانية الحرجة \(\xi_\mathrm{crit}\). ندرج قيم الجليد والسيليكات في الجدول [tab:parameters] تماشيًا مع (Kataoka2013).

نركز على سلوك الدوران بين وحدتين نتيجة هيمنة الحركة الدورانية عند عوامل الملء المنخفضة (φ<0.1) (Kataoka2013). عندما تتجاوز الإزاحة الدورانية الحد الحرج \(\xi_\mathrm{crit}\)، تتحول الحركة إلى انزياح دوار. نختار \(\xi_\mathrm{crit}=8\) Å للجليد وفقًا لـ(Kataoka2013)، وندرس تأثيره في [subsec:result:parameter]. الطاقة اللازمة للدحرجة مسافة \((\pi/2)r_0\) تعطى ب:

\[ E_\mathrm{roll} = 6\pi^2\gamma r_0\xi_\mathrm{crit}\simeq4.7\times10^{-16}\mathrm{\ J}\times\left(\frac{\gamma}{100\mathrm{\ mJ\ m^{-2}}}\right)\left(\frac{r_0}{0.1\mathrm{\ \mu m}}\right)\left(\frac{\xi_\mathrm{crit}}{8\textrm{\ \AA}}\right). \]

نضيف تخميدًا اصطناعيًا في الاتجاه العمودي بمعامل بلا بعد \(k_\mathrm{n}\) (انظر Tatsuuma2019). نختار \(k_\mathrm{n}=0.01\) وفقًا لـ(Kataoka2013)، مع العلم أن تأثير التخميد على قوة الضغط ضئيل.

إعدادات محاكاة الضغط

الخطوط العريضة لمحاكاةنا العددية كما يلي: أولًا، ننشئ تجمعًا وفق طريقة BCCA بشكل عشوائي. ثانيًا، نضغطه ببطء وبشكل متساوٍ باستخدام حدود دورية متحركة (انظر Kataoka2013 القسم 2.3).

سرعة الحدود تُعطى بـ: \[ v_\mathrm{b} = -\frac{C_\mathrm{v}}{t_\mathrm{c}}L, \] حيث \(C_\mathrm{v}\) معامل معدل الإجهاد، \(t_\mathrm{c}\) الزمن المميز (Wada2007)، و\(L\) طول الصندوق. \(t_\mathrm{c}\) تُحسب كما في النص الأصلي، ونختار \(C_\mathrm{v}=1\times10^{-7}\) كتجربة مرجعية و\(C_\mathrm{v}=3\times10^{-7}\) للمحاكاة ذات الوقت المعقول (Kataoka2013).

قياس قوة الضغط وعامل التعبئة الحجمي

نحسب قوة الضغط \(P_\mathrm{calc}\) كما في (Kataoka2013): \[ P_\mathrm{calc} = \frac{2K}{3V} + \frac{1}{3V}\left\langle\sum_{i حيث \(K\) متوسط الطاقة الحركية و\(V\) حجم الصندوق.

نعرف عامل التعبئة الحجمي بالشكل التالي: \[ \phi_\mathrm{calc} = \frac{(4/3)\pi r_0^3N}{V}. \] نأخذ المتوسط الزمني لكل 10,000 خطوة، مع اعتبار أن الخطوة الزمنية \(0.1t_\mathrm{c}\).

النتائج العددية

في هذا القسم نقدم نتائج المحاكاة لاشتقاق قوة الضغط. نجري 10 محاكاة لمجموعات أولية مختلفة ونأخذ المتوسط للتقليل من تأثير التكوِّنات العشوائية. أولًا، نعرض النتائج الأساسية في [subsec:result:fiducial]، ثم ندرس تبعيات المعاملات في [subsec:result:parameter]. (Kataoka2013) أكد أن النتائج مستقرة بالنسبة للمعاملات العددية مثل \(N\)، \(C_\mathrm{v}\)، \(k_\mathrm{n}\) (انظر الملحق [apsec:parameterdepend]).

التجارب الأساسية

تُظهر التجارب الأساسية قوة الضغط للمعاملات المرجعية. الخط المنقط يمثل الصيغة التحليلية لـ(Kataoka2013): \[ P_\mathrm{K13} = \frac{E_\mathrm{roll}}{r_0^3}\phi^3 \simeq4.7\times10^5\mathrm{\ Pa}\left(\frac{\gamma}{100\mathrm{\ mJ\ m^{-2}}}\right)\left(\frac{r_0}{0.1\mathrm{\ \mu m}}\right)^{-2}\left(\frac{\xi_\mathrm{crit}}{8\textrm{\ \AA}}\right)\phi^3. \] تتوافق المحاكاة مع هذه الصيغة عند \(\phi\lesssim0.1\)، لكن عند \(\phi>0.1\) نجد قوة ضغط أعلى بكثير وتباينًا واضحًا بين التجارب.

اعتمادات المعاملات

وجدنا أن قوة الضغط لا تعتمد على الإزاحة الدورانية الحرجة \(\xi_\mathrm{crit}\) عند \(\phi\gtrsim0.3\). أما عند \(\phi\lesssim0.3\) فتظهر تبعية قوية على \(\xi_\mathrm{crit}\) حيث تهيمن حركة الدوران. يحدث استثناء طفيف عند \(\xi_\mathrm{crit}=32\textrm{\ \AA}\) لأن حركة الالتواء تصبح مهيمنة (>16 Å).

بالنسبة للمعاملات الأخرى، تتناسب نتائج المحاكاة مع التكهن \(P\propto\gamma\,\xi_\mathrm{crit}\,r_0^{-2}\) فقط عند \(\phi\lesssim0.1\)، بينما تتباعد عند القيم الأكبر. نلاحظ تقلبات عند \(\phi<10^{-2}\) بسبب عدم ارتباط التجمعات بجميع حدود الصندوق.

المناقشة

في هذا القسم نناقش تبعيات المعاملات والفيزياء وراء قوة الضغط. أولًا، نصحح الصيغة التحليلية لعوامل تعبئة أكبر من 0.1 في [subsec:discuss:formulate]. ثانيًا، ندرس تعطيل الوحدات الأحادية في [subsec:discuss:monomerdisruption]. ثالثًا، نستعرض آليات تبديد الطاقة في [subsec:discuss:energy]. أخيرًا، نقارن نتائجنا مع الدراسات السابقة في [subsec:discuss:compare].

الصيغة المصححة لقوة الضغط

أظهرنا في القسم [sec:result] أن الصيغة البسيطة لـ(Kataoka2013) تقلل التقدير عند \(\phi>0.1\). نقترح هنا صيغة عُليا تعوّض عن حجم المونومرات المحدود:

أولًا، نعكس معادلة ([eq:Pcomp_kataoka]) لتصبح:

\[ P_\mathrm{comp} = \frac{E_\mathrm{roll}}{r_0^3}{\phi'}^3 = \frac{E_\mathrm{roll}}{r_0^3}\Bigl(\frac{NV_0}{V'}\Bigr)^3, \]

حيث \(V_0=(4/3)\pi r_0^3\) وحجم الفراغ \(V'\). ثانيًا، نطرح الحجم المستبعد \(V_\mathrm{cp}=NV_0/\phi_\mathrm{max}\) للحصول على:

\[ P_\mathrm{comp} = \frac{E_\mathrm{roll}}{r_0^3}\Bigl(\frac{NV_0}{V - NV_0/\phi_\mathrm{max}}\Bigr)^3 = \frac{E_\mathrm{roll}}{r_0^3}\Bigl(\frac{1}{\phi}-\frac{1}{\phi_\mathrm{max}}\Bigr)^{-3}, \]

مع تقريب عددي مماثل لنص المعادلة ([eq:Pcomp_mod]).

لعكس المعادلة إلى \(\phi\) كدالة لـ \(P\)، نحصل على:

\[ \phi_\mathrm{comp}=\Bigl(\frac{E_\mathrm{roll}^{1/3}}{r_0P^{1/3}}+\frac{1}{\phi_\mathrm{max}}\Bigr)^{-1}, \]

مع فرض \(\phi_\mathrm{max}=0.74\). كما نعكس معادلة ([eq:Pcomp_kataoka]) لنحصل على:

\[ \phi_\mathrm{K13}=\frac{r_0P^{1/3}}{E_\mathrm{roll}^{1/3}}. \]

تُظهر المقارنة أن الصيغة المصححة تتطابق بشكل أفضل مع نتائج المحاكاة، بخطأ أقل من 0.1 في معظم الحالات، مع ملاحظة تشوه المونومرات المرن أحيانًا.

تعطيل الوحدة الأحادية

نناقش هنا الحد الأعلى لتطبيق صيغة قوة الضغط قبل كسر الوحدات الأحادية. تشير الدراسات إلى أن الجليد ينكسر عند 5–25 ميغاباسكال عند درجات حرارة –10° إلى –20°C (Haynes1978, Petrovic2003)، بينما ينكسر زجاج السيليكا عند ~5 جيجاباسكال (Proctor1967, Bruckner1970, Kurkjian2003).

نشير إلى أن الضغط على سطح الاتصال يتركز على مساحة صغيرة \(a^2\ll r_0^2\)، حيث نصف قطر الاتصال \(a_0\) يُعطى بـ:

\[ a_0 = \Bigl(\frac{9\pi\gamma r_0^2}{4E^\ast}\Bigr)^{1/3}\simeq0.012\mathrm{\ \mu m}\Bigl(\frac{\gamma}{100\mathrm{\ mJ\ m^{-2}}}\Bigr)^{1/3}\Bigl(\frac{r_0}{0.1\mathrm{\ \mu m}}\Bigr)^{2/3}\Bigl(\frac{E^\ast}{3.7\mathrm{\ GPa}}\Bigr)^{-1/3}. \]

يمكن تقدير الحد الأعلى لقوة الضغط قبل الكسر بـ:

\[ P_\mathrm{ul}\sim P_\mathrm{dis}\Bigl(\frac{a_0}{r_0}\Bigr)^2\simeq0.014\,P_\mathrm{dis}\Bigl(\frac{\gamma}{100\mathrm{\ mJ\ m^{-2}}}\Bigr)^{2/3}\Bigl(\frac{r_0}{0.1\mathrm{\ \mu m}}\Bigr)^{-2/3}\Bigl(\frac{E^\ast}{3.7\mathrm{\ GPa}}\Bigr)^{-2/3}, \]

حيث \(P_\mathrm{dis}\) قوة التعطيل (~10 ميغاباسكال للجليد، 1 جيجاباسكال للسيليكات). بالتالي يُحسب الحد الأعلى لعامل التعبئة عبر المعادلة ([eq:phi_upperlimit]) وينتج قيمًا في النطاق 0.36–0.64 حسب المادة والحجم.

آليات تبديد الطاقة

نستعرض هنا الآليات المختلفة لتبديد الطاقة أثناء الضغط، مثل التدحرج والانزلاق والالتواء في الوصلات بين الجسيمات.

مقارنة مع الدراسات السابقة

هناك العديد من الدراسات التجريبية والعددية حول قوة الضغط شبه الثابتة لتجمعات السيليكات وثاني أكسيد السيليكون عند عوامل ملء حجمية أعلى من 0.1 (Guttler2009, Seizinger2012, Omura2017, Omura2018, Omura2021). أجرى (Seizinger2012) محاكاة عددية لضغط عمودي فقط، بينما ضغطنا كان ثلاثي الأبعاد، ورغم ذلك تتوافق النتائج العددية. أظهرت التجارب على ثاني أكسيد السيليكون (Guttler2009) اختلافًا بسبب قيمة طاقة السطح المفترضة؛ يتحسن التوافق عند افتراض \(\gamma\simeq210\pm90\mathrm{\ mJ\ m^{-2}}\). تتطلب دراسات Omura2021 لتوزيعات حجمية متعددة أيضًا \(\gamma\) أعلى للتوافق الكامل، مع بعض التباين بسبب تنظيم التصاق الجسيمات ذات الأحجام المختلفة.

ملاحظة: الخطوط المنقطة كما في اللوحة اليسرى لكن مع \(\gamma=210\mathrm{\ mJ\ m^{-2}}\).

الاستنتاجات

أجرينا محاكاة عددية لضغط تجمعات الغبار وصغنا قوة الضغط عبر كامل نطاق عوامل التعبئة الحجمية. استخدمنا نموذج التفاعل الأحادي استنادًا إلى (Dominik1997, Wada2007)، وضغطنا تجمعات BCCA ثلاثيًا الأبعاد بحدود دورية متحركة. حسبنا القوة من الطاقة الحركية ومجموع القوى على الوصلات، وأجرينا 10 محاكاة لكل مجموعة معاملات لأخذ المتوسط.

أهم النتائج:

  1. تزيد قوة الضغط شدتها بشكل حاد عند \(\phi>0.1\)، وتستقر بالنسبة للإزاحة الدورانية الحرجة عند \(\phi\gtrsim0.3\).
  2. قمنا بصياغة صيغة تحليلية مصححة تأخذ في الاعتبار التعبئة الأعلى؛ تعتمد المعادلة ([eq:Pcomp_mod]) على الفرق بين \(\phi^{-1}\) و\(\phi_\mathrm{max}^{-1}\).
  3. تهيمن حركات الالتواء والانزلاق عند التعبئة الحجمية العالية (\(\phi>0.3\))، بينما يسيطر التدحرج عند التعبئة المنخفضة (\(\phi<0.3\)).
  4. التوافق مع الدراسات السابقة جيد عدديًا، ويُحسّن التوافق التجريبي لتحميل السيليكات عند افتراض \(\gamma\simeq210\pm90\mathrm{\ mJ\ m^{-2}}\).

الدعم: منحة JSPS KAKENHI أرقام JP19J20351 وJP22J00260. استخدمنا نظام بيانات علم الفلك التابع لناسا وadstex.

اشتقاق قوة الضغط

في هذا الملحق نشرح الاشتقاق التفصيلي لقوة الضغط استنادًا إلى القسم [subsec:setting:measure]، بدءًا من معادلة الحركة ([eq:EOM_comp]) وصولًا إلى التعبير النهائي ([eq:P_compfinal]).

اعتمادات المعاملات الأخرى

في هذا الملحق نعرض عدم اعتماد قوة الضغط على عدد الوحدات الأحادية \(N\)، معامل معدل الإجهاد \(C_\mathrm{v}\)، معامل التخميد \(k_\mathrm{n}\)، وطول خطوة الزمن.