طرق لاتيس بولتزمان

تُعدّ طرق بولتزمان الشبكية (LBM) ، المُشتقة من طريقة أوتوماتا الغاز الشبكية (LGA) (نماذج هاردي- بوميو -بازيس وفريش - هاسلاشر - بوميو )، فئةً من طرق ديناميكا الموائع الحسابية (CFD) لمحاكاة الموائع . فبدلاً من حلّ معادلات نافيير-ستوكس مباشرةً، تُحاكى كثافة المائع على شبكة باستخدام عمليات التدفق والتصادم (الاسترخاء). ^{[ 1 ]} تتميز هذه الطريقة بتعدد استخداماتها ^{[ 1 ]} ، إذ يُمكن بسهولة جعل نموذج المائع يُحاكي سلوك الموائع الشائعة، مثل تعايش البخار والسائل، وبالتالي يُمكن محاكاة أنظمة الموائع مثل قطرات السائل. كما يُمكن محاكاة الموائع في بيئات معقدة، مثل الأوساط المسامية، بسهولة، بينما قد يصعب التعامل مع طرق ديناميكا الموائع الحسابية الأخرى ذات الحدود المعقدة.

على عكس طرق ديناميكا الموائع الحسابية التي تحل معادلات حفظ الخصائص الماكروية (مثل الكتلة والزخم والطاقة) عدديًا، فإن نموذج بولتزمان الشبكي (LBM) يُحاكي المائع المكون من جسيمات افتراضية، وتخضع هذه الجسيمات لعمليات انتشار وتصادم متتالية على شبكة منفصلة. وبفضل طبيعتها الجسيمية وديناميكياتها المحلية، يتميز نموذج بولتزمان الشبكي (LBM) بعدة مزايا مقارنةً بطرق ديناميكا الموائع الحسابية التقليدية الأخرى، لا سيما في التعامل مع الحدود المعقدة، ودمج التفاعلات المجهرية، وموازاة الخوارزمية. ^{[ 2 ]} ويُمكن تفسير معادلة بولتزمان الشبكية أيضًا على أنها معادلة بولتزمان للسرعة المنفصلة . وتُنتج الطرق العددية لحل نظام المعادلات التفاضلية الجزئية خريطة منفصلة، ​​يُمكن تفسيرها على أنها انتشار وتصادم الجسيمات الافتراضية.

تتضمن الخوارزمية خطوات التصادم والتدفق. تعمل هذه الخطوات على تطوير كثافة السائل.${\displaystyle \rho ({\vec {x}},t)}$، ل${\displaystyle {\vec {x}}}$المنصب و${\displaystyle t}$مع مرور الوقت. وبما أن السائل موجود على شبكة، فإن كثافته تتكون من عدد من المكونات.${\displaystyle f_{i},i=0,\ldots ,a}$يساوي عدد متجهات الشبكة المتصلة بكل نقطة من نقاط الشبكة. على سبيل المثال، تُعرض هنا متجهات الشبكة لشبكة بسيطة تُستخدم في عمليات المحاكاة ثنائية الأبعاد. يُرمز لهذه الشبكة عادةً بالرمز D2Q9، للدلالة على بُعدين وتسعة متجهات: أربعة متجهات على طول الشمال والشرق والجنوب والغرب، بالإضافة إلى أربعة متجهات إلى زوايا مربع الوحدة ، بالإضافة إلى متجه مركبتيه صفر. ثم، على سبيل المثال، المتجه${\displaystyle {\vec {e}}_{4}=(0,-1)}$أي أنها تشير إلى الجنوب تمامًا، وبالتالي ليس لها${\displaystyle x}$مكون ولكن${\displaystyle y}$مكون من${\displaystyle -1}$إذن، أحد المكونات التسعة للكثافة الكلية عند نقطة الشبكة المركزية،${\displaystyle f_{4}({\vec {x}},t)}$، هو جزء من السائل عند النقطة${\displaystyle {\vec {x}}}$يتحرك باتجاه الجنوب مباشرة، بسرعة تساوي واحدًا بوحدات الشبكة.

ثم الخطوات التي تطور السائل بمرور الوقت هي: ^{[ 1 ]}

التصادم

بالنسبة لنموذج بهاتناغار غروس وكروك (BGK) ^{[ 3 ]} ، الذي يؤدي إلى الاسترخاء إلى حالة التوازن عبر التصادمات بين جزيئات السائل، لدينا

${\displaystyle f_{i}^{\ast }({\vec {x}},t)=f_{i}({\vec {x}},t)+{\frac {f_{i}^{eq}({\vec {x}},t)-f_{i}({\vec {x}},t)}{\tau _{f}}}\,\!}$،

أين${\displaystyle f_{i}^{\ast }({\vec {x}},t)}$هي كثافة الشبكة الجديدة، و ${\displaystyle f_{i}^{eq}({\vec {x}},t)}$هي كثافة التوازن على طول الاتجاه i والتي يمكن التعبير عنها باستخدام متسلسلة تايلور (انظر أدناه، في المعادلات الرياضية للمحاكاة ):

${\displaystyle f_{i}^{eq}=\omega _{i}\rho \left(1+{\frac {3{\vec {e}}_{i}\cdot {\vec {u}}}{c^{2}}}+{\frac {9({\vec {e}}_{i}\cdot {\vec {u}})^{2}}{2c^{4}}}-{\frac {3({\vec {u}}\cdot {\vec {u}})}{2c^{2}}}\right)}$.

يفترض النموذج أن السائل يسترخي محليًا إلى حالة التوازن خلال فترة زمنية مميزة${\displaystyle \tau _{f}}$يحدد هذا المقياس الزمني اللزوجة الحركية ، فكلما كان أكبر، زادت اللزوجة الحركية.

خطوة البث

${\displaystyle f_{i}({\vec {x}}+{\vec {e}}_{i}\delta _{t},t+\delta _{t})=f_{i}^{\ast }({\vec {x}},t)\,\!}$

مثل${\displaystyle f_{i}^{\ast }({\vec {x}},t)}$هي، بحسب التعريف، كثافة السائل عند نقطة${\displaystyle {\vec {x}}}$في ذلك الوقت${\displaystyle t}$، أي يتحرك بسرعة${\displaystyle {\vec {e}}_{i}}$لكل خطوة زمنية، ثم في الخطوة الزمنية التالية${\displaystyle t+\delta _{t}}$سيكون قد تدفق إلى النقطة${\displaystyle {\vec {x}}+{\vec {e}}_{i}\delta _{t}}$.

• صُممت طريقة الشبكة البلورية (LBM) من الصفر لتعمل بكفاءة عالية على بنى متوازية ضخمة ، بدءًا من وحدات FPGA و DSP المدمجة منخفضة التكلفة وصولًا إلى وحدات معالجة الرسومات (GPUs ) والمجموعات الحاسوبية غير المتجانسة والحواسيب العملاقة (حتى مع شبكة ربط بطيئة). تُمكّن هذه الطريقة من دراسة الفيزياء المعقدة والخوارزميات المتطورة. وتؤدي الكفاءة إلى مستوى جديد نوعيًا من الفهم، إذ تسمح بحل المشكلات التي كان من المستحيل سابقًا معالجتها (أو كانت معالجتها بدقة غير كافية).
• تستند هذه الطريقة إلى وصف جزيئي للسوائل، ويمكنها دمج المصطلحات الفيزيائية مباشرةً من خلال معرفة التفاعل بين الجزيئات. ولذلك، فهي أداة لا غنى عنها في البحوث الأساسية، إذ تُقصر الدورة بين وضع النظرية وصياغة النموذج العددي المقابل.
• المعالجة المسبقة للبيانات الآلية وتوليد الشبكة في وقت يمثل جزءًا صغيرًا من إجمالي المحاكاة.
• تحليل البيانات المتوازية، والمعالجة اللاحقة، والتقييم.
• تدفق متعدد الأطوار محلول بالكامل مع قطرات وفقاعات صغيرة.
• تدفق محلول بالكامل عبر أشكال هندسية معقدة ووسائط مسامية.
• تدفق معقد ومقترن بانتقال الحرارة والتفاعلات الكيميائية.

كما هو الحال مع ديناميكا الموائع الحسابية القائمة على معادلات نافيير-ستوكس، فقد تم دمج طرق الشبكة البلورية بنجاح مع حلول خاصة بالحرارة لتمكين محاكاة انتقال الحرارة (التوصيل الحراري، والحمل الحراري، والإشعاع الحراري في المواد الصلبة). في النماذج متعددة الأطوار/المكونات، يكون سمك السطح البيني عادةً كبيرًا، ونسبة الكثافة عبره صغيرة مقارنةً بالسوائل الحقيقية. وقد تم حل هذه المشكلة مؤخرًا بواسطة يوان وشيفر، اللذين طورا نماذج شان وتشين، وسويفت، وهي، وتشين، وتشانغ. تمكنا من الوصول إلى نسب كثافة تبلغ 1000:1 بمجرد تغيير معادلة الحالة . وقد اقتُرح تطبيق تحويل غاليليو للتغلب على قيود نمذجة تدفقات الموائع عالية السرعة. ^{[ 4 ]} وقد نجحت التطورات السريعة لهذه الطريقة في محاكاة الموائع الدقيقة ، ^{[ 5 ]} ومع ذلك، لا تزال طريقة الشبكة البلورية (LBM) محدودة في محاكاة التدفقات ذات أرقام كنودسن العالية ، حيث تُستخدم طرق مونت كارلو بدلاً منها، كما أن محاكاة التدفقات ذات أرقام ماخ العالية في الديناميكا الهوائية لا تزال صعبة بالنسبة لطريقة الشبكة البلورية، ويفتقر الأمر إلى مخطط ديناميكي حراري مائي متسق. ^{[ 6 ]}

نشأت طريقة بولتزمان الشبكية (LBM) من طريقة أتمتة الغاز الشبكية (LGA)، والتي يمكن اعتبارها نموذجًا مبسطًا لديناميكيات الجزيئات الافتراضية، حيث تكون أبعاد المكان والزمان وسرعات الجسيمات منفصلة. على سبيل المثال، في نموذج FHP ثنائي الأبعاد ، ترتبط كل عقدة شبكية بجيرانها بست سرعات شبكية على شبكة مثلثية؛ ويمكن أن يوجد جسيم واحد أو لا يوجد جسيم عند كل عقدة شبكية يتحرك بسرعة شبكية معينة. بعد فترة زمنية محددة، ينتقل كل جسيم إلى العقدة المجاورة في اتجاهه؛ وتُسمى هذه العملية خطوة الانتشار أو التدفق. عندما يصل أكثر من جسيم إلى العقدة نفسها من اتجاهات مختلفة، فإنها تصطدم وتغير سرعاتها وفقًا لمجموعة من قواعد التصادم. تتناوب خطوات التدفق وخطوات التصادم. يجب أن تحافظ قواعد التصادم المناسبة على عدد الجسيمات (كتلتها) وزخمها وطاقتها قبل التصادم وبعده. تعاني خوارزميات LGA من عدة عيوب جوهرية عند استخدامها في المحاكاة الهيدروديناميكية، منها: عدم ثباتها في التدفقات السريعة، والضوضاء الإحصائية ، وضعف تناسب عدد رينولدز مع حجم الشبكة. مع ذلك، تُعدّ خوارزميات LGA مناسبة تمامًا لتبسيط نماذج تفاعل الانتشار وديناميكيات الجزيئات وتوسيع نطاقها .

كان الدافع الرئيسي للانتقال من خوارزمية LGA إلى خوارزمية LBM هو الرغبة في إزالة التشويش الإحصائي عن طريق استبدال عدد الجسيمات المنطقي في اتجاه الشبكة بمتوسطه الإحصائي، ما يُعرف بدالة توزيع الكثافة. ويصاحب هذا الاستبدال استبدال قاعدة التصادم المنفصلة بدالة متصلة تُعرف بمعامل التصادم. وفي تطوير خوارزمية LBM، يتمثل أحد التبسيطات المهمة في تقريب معامل التصادم باستخدام حد استرخاء بهاتناغار-غروس-كروك (BGK). يجعل نموذج BGK الشبكي (LBGK) عمليات المحاكاة أكثر كفاءة، ويتيح مرونة في معاملات النقل. من ناحية أخرى، ثبت أن مخطط LBM يمكن اعتباره أيضًا شكلًا منفصلًا خاصًا لمعادلة بولتزمان المتصلة. ومن خلال نظرية تشابمان-إنسكوغ ، يمكن استخلاص معادلات الاستمرارية ونافير-ستوكس الحاكمة من خوارزمية LBM.

يمكن تشغيل نماذج بولتزمان الشبكية على عدد من الشبكات المختلفة، المكعبة والمثلثة، ومع أو بدون جسيمات سكون في دالة التوزيع المنفصلة.

إحدى الطرق الشائعة لتصنيف الطرق المختلفة حسب الشبكة هي مخطط D_n Q_m . يرمز "D _n " هنا إلى " n بُعد"، بينما يرمز "Q_m " إلى " m سرعة". على سبيل المثال، D3Q15 هو نموذج بولتزمان ثلاثي الأبعاد على شبكة مكعبة، مع وجود جسيمات ثابتة. لكل عقدة شكل بلوري، ويمكنها نقل الجسيمات إلى 15 عقدة: كل من العقد الست المجاورة التي تشترك في سطح، والعقد الثماني المجاورة التي تشترك في زاوية، بالإضافة إلى العقدة نفسها. ^{[ 7 ]} (لا يحتوي نموذج D3Q15 على جسيمات تتحرك إلى العقد الـ 12 المجاورة التي تشترك في حافة؛ إضافة هذه العقد ستؤدي إلى إنشاء نموذج "D3Q27").

يجب تحويل الكميات الحقيقية، كالمكان والزمان، إلى وحدات شبكية قبل إجراء المحاكاة. أما الكميات غير البعدية، مثل عدد رينولدز ، فتبقى كما هي.

في معظم عمليات محاكاة لاتيس بولتزمان${\displaystyle \delta _{x}\,\!}$هي الوحدة الأساسية لتباعد الشبكة، لذلك إذا كان نطاق الطول${\displaystyle L\,\!}$لديه${\displaystyle N\,\!}$تُعرَّف وحدة الفضاء ببساطة على أنها وحدات شبكية على طولها بالكامل.${\displaystyle \delta _{x}=L/N\,\!}$تُعطى السرعات في محاكاة بولتزمان الشبكية عادةً بدلالة سرعة الصوت. وبالتالي، يمكن التعبير عن وحدة الزمن المنفصلة على النحو التالي:${\displaystyle \delta _{t}={\frac {\delta _{x}}{C_{s}}}\,\!}$، حيث المقام${\displaystyle C_{s}}$هي السرعة الفيزيائية للصوت. ^{[ 8 ]}

في التدفقات الصغيرة (مثل تلك التي تُرى في ميكانيكا الأوساط المسامية )، قد يؤدي استخدام سرعة الصوت الحقيقية إلى خطوات زمنية قصيرة بشكل غير مقبول. لذلك، من الشائع رفع عدد ماخ الشبكي إلى قيمة أكبر بكثير من عدد ماخ الحقيقي، والتعويض عن ذلك برفع اللزوجة أيضًا للحفاظ على عدد رينولدز . ^{[ 9 ]}

لطالما شكلت محاكاة التدفقات متعددة الأطوار/المكونات تحديًا لنظرية ديناميكا الموائع الحسابية التقليدية نظرًا لحركة وتشوه الأسطح البينية . وبشكل أساسي، تنشأ الأسطح البينية بين الأطوار المختلفة (السائل والبخار) أو المكونات (مثل الزيت والماء) من التفاعلات المحددة بين جزيئات المائع. لذا، يصعب دمج هذه التفاعلات المجهرية في معادلة نافيير-ستوكس الكلية. مع ذلك، في طريقة الشبكة البلورية (LBM)، توفر حركية الجسيمات طريقة سهلة ومتسقة نسبيًا لدمج التفاعلات المجهرية الأساسية من خلال تعديل معامل التصادم. وقد طُوّرت العديد من نماذج LBM متعددة الأطوار/المكونات. في هذه النماذج، تُولّد فواصل الأطوار تلقائيًا من ديناميكيات الجسيمات، ولا حاجة إلى معالجة خاصة للتحكم في الأسطح البينية كما في طرق ديناميكا الموائع الحسابية التقليدية. يمكن إيجاد تطبيقات ناجحة لنماذج LBM متعددة الأطوار / متعددة المكونات في أنظمة السوائل المعقدة المختلفة ، بما في ذلك عدم استقرار الواجهة، وديناميكيات الفقاعات / القطرات ، والتبلل على الأسطح الصلبة، والانزلاق البيني، والتشوهات الكهروهيدروديناميكية للقطرات.

تم مؤخراً اقتراح نموذج بولتزمان الشبكي لمحاكاة احتراق خليط الغازات، وهو نموذج قادر على استيعاب تغيرات الكثافة الكبيرة في نظام عدد ماخ المنخفض. ^{[ 10 ]}

في هذا الصدد، تجدر الإشارة إلى أنه نظرًا لأن طريقة الشبكة البلورية (LBM) تتعامل مع مجموعة أكبر من المجالات (مقارنةً بديناميكيات الموائع الحسابية التقليدية)، فإن محاكاة مخاليط الغازات التفاعلية تُطرح بعض التحديات الإضافية فيما يتعلق بمتطلبات الذاكرة، وذلك فيما يخص آليات الاحتراق التفصيلية الكبيرة. ومع ذلك، يمكن معالجة هذه المشكلات باللجوء إلى تقنيات منهجية لتقليل حجم النموذج. ^{[ 11 ] [ 12 ] [ 13 ]}

حالياً (2009)، تندرج طريقة بولتزمان الشبكية الحرارية (TLBM) ضمن إحدى الفئات الثلاث التالية: نهج السرعات المتعددة، ^{[ 14 ]} نهج الكمية القياسية السلبية، ^{[ 15 ]} وتوزيع الطاقة الحرارية. ^{[ 16 ]}

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

${\displaystyle f_{i}({\vec {x}}+{\vec {e}}_{i}\delta _{t},t+\delta _{t})=f_{i}({\vec {x}},t)+{\frac {\delta _{t}}{\tau _{f}}}(f_{i}^{eq}-f_{i}).}$

لتبسيط الأمر، اكتب${\displaystyle f_{i}({\vec {x}},t)}$مثل${\displaystyle f_{i}}$. يكون توسيع متسلسلة تايلور المبسط قليلاً كما يلي، حيث ":" هو حاصل الضرب النقطي بين الثنائيات:

${\displaystyle {\frac {\partial f_{i}}{\partial t}}+{\vec {e}}_{i}\cdot \nabla f_{i}+\left({\frac {1}{2}}{\vec {e}}_{i}{\vec {e}}_{i}:\nabla \nabla f_{i}+{\vec {e}}_{i}\cdot \nabla {\frac {\partial f_{i}}{\partial t}}+{\frac {1}{2}}{\frac {\partial ^{2}f_{i}}{\partial t^{2}}}\right)={\frac {1}{\tau }}(f_{i}^{eq}-f_{i}).}$

من خلال توسيع دالة توزيع الجسيمات إلى مكونات التوازن وعدم التوازن وباستخدام توسيع تشابمان-إنسكوج، حيث${\displaystyle K}$إذا كان عدد كنودسن هو ، فيمكن تحليل معادلة تايلور الموسعة LBE إلى مقادير مختلفة من رتبة عدد كنودسن من أجل الحصول على معادلات الاستمرارية المناسبة:

${\displaystyle f_{i}=f_{i}^{\text{eq}}+Kf_{i}^{\text{neq}},}$

${\displaystyle f_{i}^{\text{neq}}=f_{i}^{(1)}+Kf_{i}^{(2)}+O(K^{2}).}$

تخضع التوزيعات المتوازنة وغير المتوازنة للعلاقات التالية مع متغيراتها الكلية (سيتم استخدام هذه العلاقات لاحقًا، بمجرد أن تصبح توزيعات الجسيمات في "الشكل الصحيح" من أجل الانتقال من مستوى الجسيمات إلى المستوى الكلي):

${\displaystyle \rho =\sum _{i}f_{i}^{\text{eq}},}$

${\displaystyle \rho {\vec {u}}=\sum _{i}f_{i}^{\text{eq}}{\vec {e}}_{i},}$

${\displaystyle 0=\sum _{i}f_{i}^{(k)}\qquad {\text{for }}k=1,2,}$

${\displaystyle 0=\sum _{i}f_{i}^{(k)}{\vec {e}}_{i}.}$

إذن، يكون توسيع تشابمان-إنسكوج كما يلي:

${\displaystyle {\frac {\partial }{\partial t}}=K{\frac {\partial }{\partial t_{1}}}+K^{2}{\frac {\partial }{\partial t_{2}}}\qquad {\text{for }}t_{2}({\text{diffusive time-scale}})\ll t_{1}({\text{convective time-scale}}),}$

${\displaystyle {\frac {\partial }{\partial x}}=K{\frac {\partial }{\partial x_{1}}}.}$

عن طريق استبدال حالة التوازن الموسعة وحالة عدم التوازن في متسلسلة تايلور وفصلها إلى رتب مختلفة من${\displaystyle K}$، يتم اشتقاق معادلات الوسط المتصل تقريبًا.

للطلب${\displaystyle K^{0}}$:

${\displaystyle {\frac {\partial f_{i}^{\text{eq}}}{\partial t_{1}}}+{\vec {e}}_{i}\nabla _{1}f_{i}^{\text{eq}}=-{\frac {f_{i}^{(1)}}{\tau }}.}$

للطلب${\displaystyle K^{1}}$:

${\displaystyle {\frac {\partial f_{i}^{(1)}}{\partial t_{1}}}+{\frac {\partial f_{i}^{\text{eq}}}{\partial t_{2}}}+{\vec {e}}_{i}\nabla f_{i}^{(1)}+{\frac {1}{2}}{\vec {e}}_{i}{\vec {e}}_{i}:\nabla \nabla f_{i}^{\text{eq}}+{\vec {e}}_{i}\cdot \nabla {\frac {\partial f_{i}^{\text{eq}}}{\partial t_{1}}}+{\frac {1}{2}}{\frac {\partial ^{2}f_{i}^{\text{eq}}}{\partial t_{1}^{2}}}=-{\frac {f_{i}^{(2)}}{\tau }}.}$

بعد ذلك، يمكن تبسيط المعادلة الثانية باستخدام بعض العمليات الجبرية، والمعادلة الأولى إلى ما يلي:

${\displaystyle {\frac {\partial f_{i}^{\text{eq}}}{\partial t_{2}}}+\left(1-{\frac {1}{2\tau }}\right)\left[{\frac {\partial f_{i}^{(1)}}{\partial t_{1}}}+{\vec {e}}_{i}\nabla _{1}f_{i}^{(1)}\right]=-{\frac {f_{i}^{(2)}}{\tau }}.}$

بتطبيق العلاقات بين دوال توزيع الجسيمات والخصائص العيانية المذكورة أعلاه، يتم التوصل إلى معادلات الكتلة والزخم:

${\displaystyle {\frac {\partial \rho }{\partial t}}+\nabla \cdot \rho {\vec {u}}=0,}$

${\displaystyle {\frac {\partial \rho {\vec {u}}}{\partial t}}+\nabla \cdot \Pi =0.}$

موتر تدفق الزخم${\displaystyle \Pi }$ثم يكون له الشكل التالي:

${\displaystyle \Pi _{xy}=\sum _{i}{\vec {e}}_{ix}{\vec {e}}_{iy}\left[f_{i}^{eq}+\left(1-{\frac {1}{2\tau }}\right)f_{i}^{(1)}\right],}$

أين${\displaystyle {\vec {e}}_{ix}{\vec {e}}_{iy}}$هو اختصار لمربع مجموع جميع مكونات${\displaystyle {\vec {e}}_{i}}$(أي${\displaystyle \textstyle \left(\sum _{x}{\vec {e}}_{ix}\right)^{2}=\sum _{x}\sum _{y}{\vec {e}}_{ix}{\vec {e}}_{iy}}$), ويكون توزيع الجسيمات في حالة التوازن من الدرجة الثانية قابلاً للمقارنة مع معادلة نافيير-ستوكس كما يلي:

${\displaystyle f_{i}^{\text{eq}}=\omega _{i}\rho \left(1+{\frac {{\vec {e}}_{i}{\vec {u}}}{c_{s}^{2}}}+{\frac {({\vec {e}}_{i}{\vec {u}})^{2}}{2c_{s}^{4}}}-{\frac {{\vec {u}}^{2}}{2c_{s}^{2}}}\right).}$

لا يكون توزيع التوازن صالحًا إلا للسرعات المنخفضة أو أرقام ماخ المنخفضة . يؤدي إدخال توزيع التوازن مرة أخرى في موتر التدفق إلى:

${\displaystyle \Pi _{xy}^{(0)}=\sum _{i}{\vec {e}}_{ix}{\vec {e}}_{iy}f_{i}^{eq}=p\delta _{xy}+\rho u_{x}u_{y},}$

${\displaystyle \Pi _{xy}^{(1)}=\left(1-{\frac {1}{2\tau }}\right)\sum _{i}{\vec {e}}_{ix}{\vec {e}}_{iy}f_{i}^{(1)}=\nu \left(\nabla _{x}\left(\rho {\vec {u}}_{y}\right)+\nabla _{y}\left(\rho {\vec {u}}_{x}\right)\right).}$

وأخيرًا، يتم استعادة معادلة نافيير-ستوكس بافتراض أن تغير الكثافة صغير:

${\displaystyle \rho \left({\frac {\partial {\vec {u}}_{x}}{\partial t}}+\nabla _{y}\cdot {\vec {u}}_{x}{\vec {u}}_{y}\right)=-\nabla _{x}p+\nu \nabla _{y}\cdot \left(\nabla _{x}\left(\rho {\vec {u}}_{y}\right)+\nabla _{y}\left(\rho {\vec {u}}_{x}\right)\right).}$

ويستند هذا الاشتقاق إلى عمل تشين ودولين. ^{[ 17 ]}

معادلة بولتزمان المستمرة هي معادلة تطور لدالة توزيع احتمالية جسيم واحد${\displaystyle f({\vec {x}},{\vec {e}}_{i},t)}$ودالة توزيع كثافة الطاقة الداخلية${\displaystyle g({\vec {x}},{\vec {e}}_{i},t)}$(هي وآخرون) كل منهم على التوالي:

${\displaystyle \partial _{t}f+({\vec {e}}\cdot \nabla )f+F\partial _{v}f=\Omega (f),}$

${\displaystyle \partial _{t}g+({\vec {e}}\cdot \nabla )g+G\partial _{v}f=\Omega (g),}$

أين${\displaystyle g({\vec {x}},{\vec {e}}_{i},t)}$يرتبط بـ${\displaystyle f({\vec {x}},{\vec {e}}_{i},t)}$بواسطة

${\displaystyle g({\vec {x}},{\vec {e}}_{i},t)={\frac {({\vec {e}}-{\vec {u}})^{2}}{2}}f({\vec {x}},{\vec {e}}_{i},t),}$

${\displaystyle F}$هي قوة خارجية،${\displaystyle \Omega }$هو تكامل تصادمي، و${\displaystyle {\vec {e}}}$(يُشار إليه أيضًا بواسطة${\displaystyle {\vec {\xi }}}$(في الأدب) هي السرعة المجهرية. القوة الخارجية${\displaystyle F}$يرتبط ذلك بدرجة الحرارة والقوة الخارجية${\displaystyle G}$وفقًا للعلاقة أدناه. يُعد اختبار رايلي-بينارد للحمل الحراري اختبارًا نموذجيًا لنموذج المرء.${\displaystyle G}$.

${\displaystyle F={\frac {{\vec {G}}\cdot ({\vec {e}}-{\vec {u}})}{RT}}f^{\text{eq}},}$

${\displaystyle {\vec {G}}=\beta g_{0}(T-T_{avg}){\vec {k}}.}$

المتغيرات العيانية مثل الكثافة${\displaystyle \rho }$، سرعة${\displaystyle {\vec {u}}}$ودرجة الحرارة${\displaystyle T}$يمكن حسابها كعزوم لدالة توزيع الكثافة:

${\displaystyle \rho =\int f\,d{\vec {e}},}$

${\displaystyle \rho {\vec {u}}=\int {\vec {e}}f\,d{\vec {e}},}$

${\displaystyle {\frac {\rho DRT}{2}}=\rho \epsilon =\int g\,d{\vec {e}}.}$

تقوم طريقة بولتزمان الشبكية بتقسيم هذه المعادلة عن طريق حصر الفضاء في شبكة وحصر فضاء السرعة في مجموعة منفصلة من السرعات المجهرية (أي${\displaystyle {\vec {e}}_{i}=({\vec {e}}_{ix},{\vec {e}}_{iy})}$على سبيل المثال، تُعطى السرعات المجهرية في D2Q9 وD3Q15 وD3Q19 على النحو التالي:

${\displaystyle {\vec {e}}_{i}=c\times {\begin{cases}(0,0)&i=0\\(1,0),(0,1),(-1,0),(0,-1)&i=1,2,3,4\\(1,1),(-1,1),(-1,-1),(1,-1)&i=5,6,7,8\\\end{cases}}}$

${\displaystyle {\vec {e}}_{i}=c\times {\begin{cases}(0,0,0)&i=0\\(\pm 1,0,0),(0,\pm 1,0),(0,0,\pm 1)&i=1,2,...,5,6\\(\pm 1,\pm 1,\pm 1)&i=7,8,...,13,14\\\end{cases}}}$

${\displaystyle {\vec {e}}_{i}=c\times {\begin{cases}(0,0,0)&i=0\\(\pm 1,0,0),(0,\pm 1,0),(0,0,\pm 1)&i=1,2,...,5,6\\(\pm 1,\pm 1,0),(\pm 1,0,\pm 1),(0,\pm 1,\pm 1)&i=7,8,...,17,18\\\end{cases}}}$

معادلة بولتزمان المتقطعة أحادية الطور لكثافة الكتلة وكثافة الطاقة الداخلية هي:

${\displaystyle f_{i}({\vec {x}}+{\vec {e}}_{i}\delta _{t},t+\delta _{t})-f_{i}({\vec {x}},t)+F_{i}=\Omega (f),}$

${\displaystyle g_{i}({\vec {x}}+{\vec {e}}_{i}\delta _{t},t+\delta _{t})-g_{i}({\vec {x}},t)+G_{i}=\Omega (g).}$

غالبًا ما يتم تقريب عامل التصادم بواسطة عامل تصادم BGK بشرط أن يفي أيضًا بقوانين الحفظ:

${\displaystyle \Omega (f)={\frac {1}{\tau _{f}}}(f_{i}^{\text{eq}}-f_{i}),}$

${\displaystyle \Omega (g)={\frac {1}{\tau _{g}}}(g_{i}^{\text{eq}}-g_{i}).}$

في مشغل التصادم${\displaystyle f_{i}^{\text{eq}}}$هي دالة توزيع احتمالية الجسيمات المتقطعة والمتوازنة . في D2Q9 وD3Q19، تُعرض أدناه لتدفق غير قابل للانضغاط في شكله المستمر والمتقطع، حيث D و R و T هي الأبعاد وثابت الغازات العام ودرجة الحرارة المطلقة على التوالي. يُقدم الاشتقاق الجزئي لتحويل الشكل المستمر إلى شكل متقطع من خلال اشتقاق بسيط بدقة من الدرجة الثانية.

${\displaystyle f^{\text{eq}}={\frac {\rho }{(2\pi RT)^{D/2}}}e^{-{\frac {({\vec {e}}-{\vec {u}})^{2}}{2RT}}}}$

${\displaystyle ={\frac {\rho }{(2\pi RT)^{D/2}}}e^{-{\frac {({\vec {e}})^{2}}{2RT}}}e^{{\frac {{\vec {e}}{\vec {u}}}{RT}}-{\frac {{\vec {u}}^{2}}{2RT}}}}$

${\displaystyle ={\frac {\rho }{(2\pi RT)^{D/2}}}e^{-{\frac {({\vec {e}})^{2}}{2RT}}}\left(1+{\frac {{\vec {e}}{\vec {u}}}{RT}}+{\frac {({\vec {e}}{\vec {u}})^{2}}{2(RT)^{2}}}-{\frac {{\vec {u}}^{2}}{2RT}}+...\right)}$

تأجير${\displaystyle c={\sqrt {3RT}}}$ويؤدي ذلك إلى النتيجة النهائية:

${\displaystyle f_{i}^{eq}=\omega _{i}\rho \left(1+{\frac {3{\vec {e}}_{i}{\vec {u}}}{c^{2}}}+{\frac {9({\vec {e}}_{i}{\vec {u}})^{2}}{2c^{4}}}-{\frac {3({\vec {u}})^{2}}{2c^{2}}}\right)}$

${\displaystyle g^{eq}={\frac {\rho ({\vec {e}}-{\vec {u}})^{2}}{2(2\pi RT)^{D/2}}}e^{-{\frac {({\vec {e}}-{\vec {u}})^{2}}{2RT}}}}$

${\displaystyle \omega _{i}={\begin{cases}4/9&i=0\\1/9&i=1,2,3,4\\1/36&i=5,6,7,8\\\end{cases}}}$

${\displaystyle \omega _{i}={\begin{cases}1/3&i=0\\1/18&i=1,2,...,5,6\\1/36&i=7,8,...,17,18\\\end{cases}}}$

نظراً للجهود الكبيرة المبذولة في دراسة التدفق أحادي المكون، سيتم تناول نموذج TLBM التالي. كما أن نموذج TLBM متعدد المكونات/متعدد الأطوار أكثر إثارة للاهتمام وفائدة من نموذج المكون الواحد. ولمواكبة الأبحاث الحالية، يجب تحديد مجموعة جميع مكونات النظام (مثل جدران الوسط المسامي، والسوائل/الغازات المتعددة، إلخ).${\displaystyle \Psi }$مع العناصر${\displaystyle \sigma _{j}}$.

${\displaystyle f_{i}^{\sigma }({\vec {x}}+{\vec {e}}_{i}\delta _{t},t+\delta _{t})-f_{i}^{\sigma }({\vec {x}},t)+F_{i}={\frac {1}{\tau _{f}^{\sigma }}}(f_{i}^{\sigma ,eq}(\rho ^{\sigma },v^{\sigma })-f_{i}^{\sigma })}$

معامل الاسترخاء،${\displaystyle \tau _{f}^{\sigma _{j}}\,\!}$، يرتبط باللزوجة الحركية ،${\displaystyle \nu _{f}^{\sigma _{j}}\,\!}$، وفقًا للعلاقة التالية:

${\displaystyle \nu _{f}^{\sigma _{j}}=(\tau _{f}^{\sigma _{j}}-0.5)c_{s}^{2}\delta _{t}.}$

لحظات​​${\displaystyle f_{i}\,\!}$أوجد الكميات المحفوظة المحلية. الكثافة معطاة بواسطة

${\displaystyle \rho =\sum _{\sigma }\sum _{i}f_{i}\,\!}$

${\displaystyle \rho \epsilon =\sum _{i}g_{i}\,\!}$

${\displaystyle \rho ^{\sigma }=\sum _{i}f_{i}^{\sigma }\,\!}$

ومتوسط ​​السرعة المرجح،${\displaystyle {\vec {u'}}\,\!}$ويتم إعطاء الزخم المحلي بواسطة

${\displaystyle {\vec {u'}}=\left(\sum _{\sigma }{\frac {\rho ^{\sigma }{\vec {u^{\sigma }}}}{\tau _{f}^{\sigma }}}\right)/\left(\sum _{\sigma }{\frac {\rho ^{\sigma }}{\tau _{f}^{\sigma }}}\right)}$

${\displaystyle \rho ^{\sigma }{\vec {u^{\sigma }}}=\sum _{i}f_{i}^{\sigma }{\vec {e}}_{i}.}$

${\displaystyle v^{\sigma }={\vec {u'}}+{\frac {\tau _{f}^{\sigma }}{\rho ^{\sigma }}}{\vec {F}}^{\sigma }}$

في المعادلة أعلاه لسرعة التوازن${\displaystyle v^{\sigma }\,\!}$، ال${\displaystyle {\vec {F}}^{\sigma }\,\!}$يمثل هذا المصطلح قوة التفاعل بين أحد المكونات والمكونات الأخرى. ولا يزال موضع نقاش واسع النطاق، إذ يُعد عادةً مُعامل ضبط يُحدد كيفية تفاعل السوائل مع بعضها البعض، أو مع الغازات، وما إلى ذلك. وقد أورد فرانك وآخرون نماذج حالية لهذا المصطلح. ومن بين الاشتقاقات الشائعة الاستخدام: نموذج غونستنسن الكروموديناميكي، ومنهج سويفت القائم على الطاقة الحرة لأنظمة السائل/البخار والسوائل الثنائية، ونموذج هي القائم على التفاعل بين الجزيئات، ومنهج إينامورو، ومنهج لي ولين. ^{[ 18 ]}

فيما يلي الوصف العام لـ${\displaystyle {\vec {F}}^{\sigma }\,\!}$كما ورد في العديد من المؤلفين. ^{[ 19 ] [ 20 ]}

${\displaystyle {\vec {F}}^{\sigma }=-\psi ^{\sigma }({\vec {x}})\sum _{\sigma _{j}}H^{\sigma \sigma _{j}}({\vec {x}},{\vec {x}}')\sum _{i}\psi ^{\sigma _{j}}({\vec {x}}+{\vec {e}}_{i}){\vec {e}}_{i}\,\!}$

${\displaystyle \psi ({\vec {x}})\,\!}$الكتلة الفعالة و${\displaystyle H({\vec {x}},{\vec {x}}')\,\!}$دالة غرين تمثل التفاعل بين الجسيمات مع${\displaystyle {\vec {x}}'\,\!}$مثل الموقع المجاور. مُرضٍ${\displaystyle H({\vec {x}},{\vec {x}}')=H({\vec {x}}',{\vec {x}})\,\!}$وأين${\displaystyle H({\vec {x}},{\vec {x}}')>0\,\!}$يمثل هذا قوى تنافر. بالنسبة لـ D2Q9 و D3Q19، يؤدي هذا إلى

${\displaystyle H^{\sigma \sigma _{j}}({\vec {x}},{\vec {x}}')={\begin{cases}h^{\sigma \sigma _{j}}&\left|{\vec {x}}-{\vec {x}}'\right|\leq c\\0&\left|{\vec {x}}-{\vec {x}}'\right|>c\\\end{cases}}}$

${\displaystyle H^{\sigma \sigma _{j}}({\vec {x}},{\vec {x}}')={\begin{cases}h^{\sigma \sigma _{j}}&\left|{\vec {x}}-{\vec {x}}'\right|=c\\h^{\sigma \sigma _{j}}/2&\left|{\vec {x}}-{\vec {x}}'\right|={\sqrt {2c}}\\0&{\text{otherwise}}\\\end{cases}}}$

تستخدم الكتلة الفعالة، كما اقترحها شان وتشن، الكتلة الفعالة التالية لنظام أحادي المكون ومتعدد الأطوار . كما تُعطى معادلة الحالة في حالة النظام أحادي المكون ومتعدد الأطوار.

${\displaystyle \psi ({\vec {x}})=\psi (\rho ^{\sigma })=\rho _{0}^{\sigma }\left[1-e^{(-\rho ^{\sigma }/\rho _{0}^{\sigma })}\right]\,\!}$

${\displaystyle p=c_{s}^{2}\rho +c_{0}h[\psi ({\vec {x}})]^{2}\,\!}$

حتى الآن، يبدو أن${\displaystyle \rho _{0}^{\sigma }\,\!}$و${\displaystyle h^{\sigma \sigma _{j}}\,\!}$هي ثوابت حرة يمكن ضبطها، ولكن بمجرد إدخالها في معادلة حالة النظام ، يجب أن تحقق العلاقات الديناميكية الحرارية عند النقطة الحرجة بحيث${\displaystyle (\partial P/\partial {\rho })_{T}=(\partial ^{2}P/\partial {\rho ^{2}})_{T}=0\,\!}$و${\displaystyle p=p_{c}\,\!}$بالنسبة لنظام التشغيل EOS،${\displaystyle c_{0}\,\!}$تبلغ قيمتها 3.0 بالنسبة لـ D2Q9 و D3Q19 بينما تساوي 10.0 بالنسبة لـ D3Q15. ^{[ 21 ]}

أظهر يوان وشيفر لاحقًا ^{[ 22 ]} أن كثافة الكتلة الفعالة تحتاج إلى تغيير لمحاكاة التدفق متعدد الأطوار بدقة أكبر. قارنا معادلات الحالة شان وتشن (SC)، وكارناهان-ستارلينغ (C–S)، وفان دير فالس (vdW)، وريدليش-كوانغ (R–K)، وريدليش-كوانغ سواف (RKS)، وبينغ-روبنسون (P–R). وكشفت نتائجهم أن معادلة الحالة SC غير كافية، وأن معادلات الحالة C–S، وP–R، وR–K، وRKS أكثر دقة في نمذجة التدفق متعدد الأطوار لمكون واحد.

بالنسبة لطرق لاتيس بولتزمان متساوية الحرارة الشائعة، تُعد هذه الكميات هي الكميات المحفوظة الوحيدة. أما النماذج الحرارية فتحافظ أيضاً على الطاقة، وبالتالي تحتوي على كمية محفوظة إضافية.

${\displaystyle \rho \theta +\rho uu=\sum _{i}f_{i}{\vec {e}}_{i}{\vec {e}}_{i}.}$

عادةً ما يتم تطبيق طرق بولتزمان الشبكية على الشبكات المنتظمة، ومع ذلك فإن استخدام الشبكة غير المنظمة يمكن أن يساعد في حل الحدود المعقدة، وتتكون الشبكات غير المنظمة من مثلثات أو رباعيات الأوجه مع اختلافات.

بافتراض${\displaystyle \Omega ^{j}}$هو حجم يتكون من جميع مراكز ثقل رباعيات الأوجه، والوجوه والحواف المتصلة بالرأس${\displaystyle {\boldsymbol {v}}^{j}}$دالة كثافة السرعة المنفصلة:

${\displaystyle f_{i}({\boldsymbol {v}}^{j},t+\delta t)=f_{i}({\boldsymbol {v}}^{j},t)-\delta t\sum _{k}S_{i}^{jk}f_{i}({\boldsymbol {v}}^{k},t)-{\delta t \over \tau }\sum _{k}C^{jk}(f_{i}({\boldsymbol {v}}^{k},t)-f_{i}^{eq}({\boldsymbol {v}}^{k}))}$

أين${\displaystyle {\boldsymbol {v}}^{k}}$هي موضع الرأس وجيرانه، و:

${\displaystyle C^{jk}={1 \over V^{j}}\int _{\Omega ^{j}}w_{k}({\boldsymbol {x}})d\Omega }$

${\displaystyle S_{i}^{jk}={1 \over V^{j}}\oint _{\partial \Omega ^{j}}({\vec {e_{i}}}{\vec {n}})w_{k}({\boldsymbol {x}})d\Omega }$

أين${\displaystyle w_{k}({\boldsymbol {x}})}$هي أوزان الاستيفاء الخطي لـ${\displaystyle {\boldsymbol {x}}}$بواسطة رؤوس المثلث أو رباعيات الأوجه التي${\displaystyle {\boldsymbol {x}}}$^{[ 23 ]}​

خلال السنوات الأخيرة، أثبتت طريقة الشبكة البلورية (LBM) أنها أداة فعّالة لحل المشكلات على مختلف المقاييس الزمنية والطولية. ومن تطبيقاتها:

• تدفقات الوسائط المسامية ^{[ 24 ]}
• التدفقات الطبية الحيوية
• علوم الأرض (ترشيح التربة).
• علوم الطاقة (خلايا الوقود ^{[ 25 ]} ).

هذا تطبيق أساسي لـ LBM على شبكة 100x100، باستخدام لغة بايثون :

# هذا برنامج محاكاة للسوائل يستخدم طريقة لاتيس بولتزمان. # يستخدم D2Q9 والحدود الدورية، ولم يستخدم أي مكتبة خارجية. # يُولّد تموجين عند 50،50 و50،40. # المرجع: رسالة الماجستير لإرلند ماغنوس فيجن، "طريقة لاتيس بولتزمان مع تطبيقات في الصوتيات". # لويكيبيديا بموجب ترخيص CC-BY-SA. استيراد الرياضيات# تعريف بعض الدوال المساعدة def sum ( a ): s = 0 for e in a : s = s + e return s# الأوزان في D2Q9 الأوزان = [ 1 / 36 , 1 / 9 , 1 / 36 , 1 / 9 , 4 / 9 , 1 / 9 , 1 / 36 , 1 / 9 , 1 / 36 ] # متجهات السرعة المنفصلة متجهات السرعة المنفصلة = [ [ - 1 , 1 ], [ 0 , 1 ], [ 1 , 1 ], [ - 1 , 0 ], [ 0 , 0 ], [ 1 , 0 ], [ - 1 , - 1 ], [ 0 , - 1 ], [ 1 , - 1 ], ]# فئة Field2D class Field2D : def __init__ ( self , res : int ) : self . field = [] for b in range ( res ): fm = [] for a in range ( res ): fm . append ([ 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 ]) self . field . append ( fm [:]) self . res = res# هذا يعرض المحاكاة، ولا يمكن استخدامه إلا في طرفية @staticmethod def VisualizeField ( a , sc , res ): stringr = "" for u in range ( res ): row = "" for v in range ( res ): n = int ( u * a . res / res ) x = int ( v * a . res / res ) vx = velocityField [ n ][ x ][ 0 ] vy = velocityField [ n ][ x ][ 1 ] r = max ( 0 , min ( 255 , int ( 127 + sc * vx ))) g = max ( 0 , min ( 255 , int ( 127 + sc * vy ))) col = " \033 [38;2; {0} ; {1} ; {2} m██" . تنسيق ( r ، g ، 0 ) صف = صف + عمود طباعة ( صف ) سلسلة نصية = سلسلة نصية + صف + " \n " إرجاع سلسلة نصية# زخم المجال def Momentum ( self , x , y ): return velocityField [ y ][ x ][ 0 ] * sum ( self . field [ y ][ x ]), velocityField [ y ][ x ][ 1 ] * sum ( self . field [ y ][ x ])# دقة المحاكاة res = 100 a = Field2D ( res ) # حقل السرعة velocityField = [] for DummyVariable in range ( res ): DummyList = [] for DummyVariable2 in range ( res ): DummyList . append ([ 0 , 0 ]) velocityField . append ( DummyList [:]) # حقل الكثافة DensityField = [] for DummyVariable in range ( res ): DummyList = [] for DummyVariable2 in range ( res ): DummyList . append ( 1 ) DensityField . أضف ( قائمة وهمية [:]) # حدد الشرط الابتدائي DensityField [ 50 ][ 50 ] = 2 DensityField [ 40 ][ 50 ] = 2 # الحد الأقصى لخطوات الحل MaxSteps = 120 # سرعة الصوت، تحديدًا 1/√3 ≈ 0.57 SpeedOfSound = 1 / math.sqrt ( 3 ) # ثابت استرخاء الزمن TimeRelaxationConstant = 0.5 # حل for s in range ( MaxSteps ) : # خطوة التصادم df = Field2D ( res ) for y in range ( res ): for x in range ( res ): for v in range ( 9 ): Velocity = a.field [ y ] [ x ][ v ] FirstTerm = Velocity # سرعة التدفق FlowVelocity = velocityField[ y ][ x ] Dotted = ( FlowVelocity [ 0 ] * DiscreteVelocityVectors [ v ][ 0 ] + FlowVelocity [ 1 ] * DiscreteVelocityVectors [ v ][ 1 ] ) # # تفسير تايلور لحد التوازن taylor = ( 1 + (( Dotted ) / ( SpeedOfSound ** 2 )) + (( Dotted ** 2 ) / ( 2 * SpeedOfSound ** 4 )) - ( ( FlowVelocity [ 0 ] ** 2 + FlowVelocity [ 1 ] ** 2 ) / ( 2 * SpeedOfSound ** 2 ) ) ) # كثافة التيار density = DensityField [ y ][ x ] # التوازن equilibrium = density * taylor * Weights [ v ] SecondTerm = ( equilibrium - السرعة ) / ثابت استرخاء الزمن df.field [ y ][ x ][ v ] = الحد الأول + الحد الثاني # خطوة البث for y in range ( 0 , res ) : for x in range ( 0 , res ): for v in range ( 9 ): # الهدف، نقطة الشبكة التي تحلها هذه التكرارات TargetY = y + DiscreteVelocityVectors [ v ][ 1 ] TargetX = x + DiscreteVelocityVectors [ v ][ 0 ]# حدود دورية إذا كان TargetY == res و TargetX == res : a.field [ TargetY - res ] [ TargetX - res ][ v ] = df.field [ y ] [ x ][ v ] وإذا كان TargetX == res : a.field [ TargetY ] [ TargetX - res ] [ v ] = df.field [ y ] [ x ] [ v ] وإذا كان TargetY == res : a.field [ TargetY - res ] [ TargetX ][ v ] = df.field [ y ] [ x ] [ v ] وإذا كان TargetY == -1 و TargetX == -1 : a.field [ TargetY + res ] [ TargetX + res ] [ v ] = df.field [ y ] [ x ] [ v ] وإذا كان TargetX == -1 : a . إذا كان TargetY يساوي -1 ، فإن الحقل [ TargetY + res ] [ TargetX + v ] يساوي الحقل [ y ] [ x ] [ v ] في إطار البيانات . وإلا ، فإن الحقل [ TargetY ] يساوي الحقل [ y ] [ x ] [ v ] في إطار البيانات .​​​TargetX ][ v ] = df . field [ y ] [ x ][ v ] # حساب المتغيرات الكلية لـ y في النطاق ( res ): لـ x في النطاق ( res ): # إعادة حساب حقل الكثافة DensityField [ y ][ x ] = sum ( a.field [ y ][ x ]) # إعادة حساب سرعة التدفق FlowVelocity = [ 0 , 0 ] لـ DummyVariable في النطاق ( 9 ): FlowVelocity [ 0 ] = ( FlowVelocity [ 0 ] + DiscreteVelocityVectors [ DummyVariable ][ 0 ] * a.field [ y ] [ x ] [ DummyVariable ] ) لـ DummyVariable في النطاق ( 9 ): FlowVelocity [ 1 ] = ( FlowVelocity [ 1 ] + DiscreteVelocityVectors [ DummyVariable ] [ 1 ] * a.field [ y ][ x ][ DummyVariable ] ) FlowVelocity [ 0 ] = FlowVelocity [ 0 ] / DensityField [ y ][ x ] FlowVelocity [ 1 ] = FlowVelocity [ 1 ] / DensityField [ y ] [ x ] # إدراج في حقل السرعة velocityField [ y ][ x ] = FlowVelocity # عرض الحقل ثنائي الأبعاد . VisualizeField ( a , 5000 , 100))