انتقل إلى المحتوى

طريقة Newton-Raphson

الخطية وعلاقة التكرار

يعطي الشغل الافتراضي للقوى الخارجية وتجميع المعادلة العامة معادلة غير خطية للإزاحة العقدية عند الزمن \(t_{n+1}\)، أي \(\boldsymbol{u}_{n+1}\)، وتُحل بطريقة Newton-Raphson. ويُفترض أن الإزاحة العقدية حتى الزمن \(t_n\)، أي \(\boldsymbol{u}_n\)، معلومة، ويُتخذ تزايد الإزاحة \(\Delta\boldsymbol{u}\) متغيرًا مجهولًا يُراد تحديده

\[ \boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u} \]

فيما يلي، يُهمَل اعتماد متجه القوة الخارجية على الإزاحة العقدية، ومع \(\boldsymbol{F}(\boldsymbol{u}_{n+1}) = \boldsymbol{F}_{n+1}\)،

\[ \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u}) - \boldsymbol{F}_{n+1} = \boldsymbol{0} \]

تُحل المعادلة.

عند الحل الحالي \(\Delta\boldsymbol{u}\)، نعرّف الصلابة المماسية

\[ \boldsymbol{K} = \left. \frac{\partial \boldsymbol{Q}}{\partial \boldsymbol{u}} \right|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}} \]

وباستخدامها، تعطي خطية المعادلة غير الخطية

\[ \boldsymbol{K}\, d\boldsymbol{u} + \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u}) - \boldsymbol{F}_{n+1} = \boldsymbol{0} \]

لنفترض أن مقدار التصحيح في التكرار رقم \(i\) هو \(d\boldsymbol{u}_i\)، وأن متجه الباقي في بداية التكرار هو

\[ \boldsymbol{R}_{i-1} = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u}) \]

وعندئذ تكون علاقة التكرار

\[ \boldsymbol{K}_i\, d\boldsymbol{u}_i = \boldsymbol{R}_{i-1}, \qquad \Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i \]

ومن ثم فإن الباقي \(\boldsymbol{R}_i\) كمية تقابل عدم اتزان القوى عن حالة الاتزان.

تكوين مصفوفة الصلابة المماسية

تُنشأ الصلابة المماسية \(\boldsymbol{K} = \partial\boldsymbol{Q}/\partial\boldsymbol{u}\) بأخذ المشتقة الجزئية لمتجه القوة الداخلية للعنصر، المحسوب في تقسيم الشغل الافتراضي للقوى الداخلية، بالنسبة إلى الإزاحة العقدية، ثم تكامل دوال التكامل الناتجة على مستوى العنصر فوق مجال كل عنصر وتجميعها. وإذا رمزنا إلى دالة التكامل على مستوى العنصر بـ \(\boldsymbol{K}^e_X\) (ترميز الوضع المرجعي، صياغة TL) أو \(\boldsymbol{K}^e_x\) (ترميز الوضع الحالي، صياغة UL)، فإن الصلابة المماسية للعنصر هي

\[ \boldsymbol{K}^e = \int_{\Omega^e_0} \boldsymbol{K}^e_X\, dV \quad (\text{TL}), \qquad \boldsymbol{K}^e = \int_{\Omega^e} \boldsymbol{K}^e_x\, dv \quad (\text{UL}) \]

فيما يلي الصيغ النهائية لدوال تكامل TL/UL. وفي كلتا الحالتين تُحلل إلى مجموع حد الصلابة المادية (حد الإزاحة الابتدائية) وحد الصلابة الهندسية (حد الإجهاد الابتدائي).

صياغة Total Lagrange

في صياغة Total Lagrange، يُفترض وجود علاقة خطية بين معدل إجهاد Piola-Kirchhoff الثاني \(\dot{\boldsymbol{S}}\) ومعدل انفعال Green-Lagrange \(\dot{\boldsymbol{E}}\)، أي \(\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}\). وهذا يقابل قوانين السلوك للمواد المرنة الخطية (مواد St. Venant-Kirchhoff) والمواد فائقة المرونة، ويستخدم FrontISTR صياغة Total Lagrange لهذه المواد. وعندئذ تُكتب دالة تكامل الصلابة المماسية للعنصر بصيغة موترية كما يلي

\[ \delta\boldsymbol{u}^{eT}\, \boldsymbol{K}^e_X\, \dot{\boldsymbol{u}}^e = \dot{\boldsymbol{S}}:\delta\boldsymbol{E} + \boldsymbol{S}:(\delta\boldsymbol{F}^T \dot{\boldsymbol{F}}) \]

الحد الأول في الطرف الأيمن هو حد الصلابة المادية (حد الإزاحة الابتدائية)، والحد الثاني هو حد الصلابة الهندسية (حد الإجهاد الابتدائي).

في تنفيذ FrontISTR، تُقيَّم دالة التكامل هذه بصيغة مصفوفية باستخدام ترميز Voigt:

\[ \boldsymbol{K}^e_X = (\boldsymbol{B}_L + \boldsymbol{B}_{NL})^T\, \tilde{\boldsymbol{C}}\, (\boldsymbol{B}_L + \boldsymbol{B}_{NL}) + \boldsymbol{F}_9^T\, \boldsymbol{S}_9\, \boldsymbol{F}_9 \]

المصفوفات كما يلي. \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}\) هما مصفوفتا B المقدمتان في تقسيم الشغل الافتراضي للقوى الداخلية، و\(\tilde{\boldsymbol{C}}\) هو تمثيل Voigt لموتر قانون السلوك \(\boldsymbol{\mathsf{C}}\)، أي مصفوفة صلابة مادية \(6\times 6\) (ترميز الموترات والأسس الرياضية). أما \(\boldsymbol{S}_9, \boldsymbol{F}_9\) فهما مصفوفتا إعادة الترتيب التاليتان المستخدمتان للتعبير عن حد الصلابة الهندسية في صورة حاصل ضرب مصفوفات. أولًا، لموتر من الرتبة الثانية \(\boldsymbol{A}\) بحجم \(3\times 3\)، نعرّف الترميز \([\,\cdot\,]\) الذي يعيد ترتيبه في متجه ذي 9 مركبات كما يلي

\[ [\boldsymbol{A}] = (A_{11}, A_{21}, A_{31}, A_{12}, A_{22}, A_{32}, A_{13}, A_{23}, A_{33})^T \]

وفق هذا التعريف، تعبّر \(\boldsymbol{F}_9\) عن تغير تدرج التشوه بالصورة \([\delta\boldsymbol{F}] = \boldsymbol{F}_9\, \delta\boldsymbol{u}^e\)، وهي مصفوفة \(9\times d n_e\). ولكل عقدة عنصر \(\alpha = 1, \ldots, n_e\)، تكون الكتلة المقابلة \(9\times d\) هي

\[ [\boldsymbol{F}_9]_\alpha = \begin{bmatrix} (\partial N_\alpha^e/\partial X_1)\, \boldsymbol{I} \\ (\partial N_\alpha^e/\partial X_2)\, \boldsymbol{I} \\ (\partial N_\alpha^e/\partial X_3)\, \boldsymbol{I} \end{bmatrix} \qquad (\boldsymbol{I} \text{ هي } 3\times 3 \text{ مصفوفة الوحدة}) \]

وتُعطى \(\boldsymbol{F}_9 = [[\boldsymbol{F}_9]_1, \ldots, [\boldsymbol{F}_9]_{n_e}]\)، مع ترتيب الكتل أفقيًا وفق ترتيب عقد العنصر. وتُختار \(\boldsymbol{S}_9\) بحيث يُعبَّر، عند جمعها مع هذه المصفوفة، عن حد الصلابة الهندسية بالشكل \(\delta\boldsymbol{u}^{eT}\, \boldsymbol{F}_9^T \boldsymbol{S}_9 \boldsymbol{F}_9\, \dot{\boldsymbol{u}}^e\)؛ وهي المصفوفة التالية \(9\times 9\)

\[ \boldsymbol{S}_9 = \begin{bmatrix} S_{11} \boldsymbol{I} & S_{12} \boldsymbol{I} & S_{13} \boldsymbol{I} \\ S_{21} \boldsymbol{I} & S_{22} \boldsymbol{I} & S_{23} \boldsymbol{I} \\ S_{31} \boldsymbol{I} & S_{32} \boldsymbol{I} & S_{33} \boldsymbol{I} \end{bmatrix} \]

وهذه هي المصفوفة الناتجة.

صياغة Updated Lagrange

في صياغة Updated Lagrange، يُفترض وجود علاقة خطية بين معدل Jaumann لموتر إجهاد Kirchhoff النسبي \(\hat{\boldsymbol{\sigma}}^{\nabla J}\) وموتر معدل التشوه \(\boldsymbol{D}\)، أي \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\). وهذه هي صيغة قانون سلوك hypoelastic شائعة للمواد المرنة الخطية واللدنة المرنة ومواد الزحف، ويستخدم FrontISTR صياغة Updated Lagrange لهذه المواد. وعندئذ تُكتب دالة تكامل الصلابة المماسية للعنصر، معبرًا عنها في الوضع الحالي، بصيغة موترية كما يلي

\[ \delta\boldsymbol{u}^{eT}\, \boldsymbol{K}^e_x\, \dot{\boldsymbol{u}}^e = \boldsymbol{\sigma}^{\nabla T}:\delta\boldsymbol{A}_{(L)} + \boldsymbol{\sigma}:(\delta\boldsymbol{F}_t^T\, \boldsymbol{L}) \]

حيث إن \(\boldsymbol{\sigma}^{\nabla T}\) هو معدل Truesdell، و\(\boldsymbol{A}_{(L)}\) هو الجزء الخطي من انفعال Almansi، و\(\boldsymbol{F}_t = \partial\boldsymbol{u}/\partial\boldsymbol{x}\) هو تدرج الإزاحة بالنسبة إلى الوضع الحالي، و\(\boldsymbol{L}\) هو موتر تدرج السرعة. الحد الأول في الطرف الأيمن هو حد الصلابة المادية، والحد الثاني هو حد الصلابة الهندسية.

في تنفيذ FrontISTR، تُقيَّم دالة التكامل هذه بصيغة مصفوفية باستخدام ترميز Voigt:

\[ \boldsymbol{K}^e_x = \boldsymbol{b}^T\, (\tilde{\boldsymbol{C}} - \boldsymbol{G})\, \boldsymbol{b} + \boldsymbol{f}_9^T\, \boldsymbol{\sigma}_9\, \boldsymbol{f}_9 \]

هنا، \(\boldsymbol{b}\) هي مصفوفة B المكوَّنة في الوضع الحالي (تقسيم الشغل الافتراضي للقوى الداخلية). أما \(\boldsymbol{\sigma}_9, \boldsymbol{f}_9\) فتُستخرجان من \(\boldsymbol{S}_9, \boldsymbol{F}_9\) المعرّفتين لصياغة TL باستبدال إجهاد PK الثاني \(\boldsymbol{S}\) بإجهاد Cauchy \(\boldsymbol{\sigma}\)، وتدرج الوضع المرجعي \(\partial N_\alpha^e/\partial X_i\) بتدرج الوضع الحالي \(\partial N_\alpha^e/\partial x_i\).

\(\boldsymbol{G}\) هي مصفوفة تصحيح تعتمد على إجهاد Cauchy، وتلزم لجعل قانون السلوك hypoelastic \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) متسقًا مع إطار الصلابة المماسية بوصفه قانون سلوك قائمًا على معدل Truesdell. وتُستخرج بترتيب مركبات الموتر من الرتبة الرابعة \(G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}\) في تمثيل Voigt بحجم \(6\times 6\) كما يلي

\[ \boldsymbol{G} = \begin{bmatrix} 2\sigma_{11} & 0 & 0 & \sigma_{12} & 0 & \sigma_{31} \\ 0 & 2\sigma_{22} & 0 & \sigma_{12} & \sigma_{23} & 0 \\ 0 & 0 & 2\sigma_{33} & 0 & \sigma_{23} & \sigma_{31} \\ \sigma_{12} & \sigma_{12} & 0 & \tfrac{\sigma_{11}+\sigma_{22}}{2} & \tfrac{\sigma_{12}}{2} & \tfrac{\sigma_{23}}{2} \\ 0 & \sigma_{23} & \sigma_{23} & \tfrac{\sigma_{12}}{2} & \tfrac{\sigma_{22}+\sigma_{33}}{2} & \tfrac{\sigma_{12}}{2} \\ \sigma_{31} & 0 & \sigma_{31} & \tfrac{\sigma_{23}}{2} & \tfrac{\sigma_{12}}{2} & \tfrac{\sigma_{33}+\sigma_{11}}{2} \end{bmatrix} \]

وهذه هي المصفوفة الناتجة.

تجميع مصفوفة الصلابة العامة

تُستخرج الصلابة المماسية العامة \(\boldsymbol{K}\) بتقسيم صلابة كل عنصر \(\boldsymbol{K}^e\) إلى كتل \(d\times d\) هي \(\boldsymbol{K}^e_{\alpha\beta}\) لكل زوج من العقد، وباستخدام مجموعة التجميع ذات الموتر من الرتبة الثانية \(\mathcal{E}^2(i_g, i_h)\) المقدمة في تجميع الكميات الفيزيائية لعقد العناصر:

\[ \boldsymbol{K}_{i_gi_h} = \sum_{(e,\alpha,\beta) \in \mathcal{E}^2(i_g, i_h)} \boldsymbol{K}^e_{\alpha\beta} \]

تُرتب القيم الناتجة في مصفوفة يكون صفها \(i_g\) وعمودها \(i_h\). وفي التنفيذ، لا تُنشأ المجموعة \(\mathcal{E}^2\) صراحة؛ بل تُضاف الكتل المقابلة مباشرة داخل حلقة العناصر. والمصفوفة مربعة، وبُعدها يساوي درجات الحرية لكل عقدة \(\times\) العدد الكلي للعقد \(n_g\)، ولكن لأن المركبات بين العقد غير المتصلة عبر عناصر تساوي \(0\)، فإنها تُخزن بصيغة مصفوفة متناثرة.

لمصفوفات صلابة العناصر في صيغتي TL وUL الشكل نفسه، باستثناء تبديل الوضع المرجعي (إحداثيات العقد والمصدر المستخدم لتكوين مصفوفة B) ووجود مصفوفة \(\boldsymbol{G}\) أو عدمه. لذلك ينفذ FrontISTR الصيغتين في روتين فرعي مشترك.

خوارزمية التكرار

تلخيصًا لما سبق، في بداية التكرار نضع \(\Delta\boldsymbol{u} = \boldsymbol{0}\) ونحسب الباقي الابتدائي \(\boldsymbol{R}_0 = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n)\). ثم، في التكرار رقم \(i\)، تُنفذ الخطوات التالية.

  1. عند الإزاحة الحالية \(\boldsymbol{u}_n + \Delta\boldsymbol{u}\)، احسب الصلابة المماسية \(\boldsymbol{K}_i = \partial\boldsymbol{Q}/\partial\boldsymbol{u}|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}}\) باستخدام الإجراء الموضح في تكوين مصفوفة الصلابة المماسية.
  2. لتطبيق الشروط الحدية الهندسية، عدّل مصفوفة الصلابة المماسية ومتجه الباقي لدرجات الحرية الخاضعة لقيود الإزاحة، للحصول على \(\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}\) (معالجة الشروط الحدية الهندسية).
  3. حل المعادلة الخطية \(\tilde{\boldsymbol{K}}_i\, d\boldsymbol{u}_i = \tilde{\boldsymbol{R}}_{i-1}\) للحصول على مقدار التصحيح \(d\boldsymbol{u}_i\). وغالبًا ما تمثل هذه العملية معظم الكلفة الحسابية للحساب التكراري.
  4. حدّث تزايد الإزاحة إلى \(\Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i\)، واحسب تبعًا لذلك متجه القوة الداخلية \(\boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\) والباقي \(\boldsymbol{R}_i = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\).
  5. افحص التقارب وأنهِ التكرار إذا تحقق. تظهر في الباقي \(\boldsymbol{R}_i\)، عند درجات الحرية الخاضعة للشروط الحدية الهندسية، مركبات تقابل ردود أفعال القيود؛ لذلك يُنشأ مؤشر التقارب من \(\tilde{\boldsymbol{R}}_i\) بعد استبعاد هذه المركبات. وتُعرض مؤشرات التقارب والعتبات المحددة في معايير التقارب. وإذا لم يتحقق التقارب وبلغ الحساب حد التكرار، يُعد التكرار فاشلًا.

عند تقارب التكرار، أضف \(\Delta\boldsymbol{u}\) المتقاربة إلى \(\boldsymbol{u}_n\) للحصول على الإزاحة التراكمية عند الزمن \(t_{n+1}\)، أي \(\boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u}\)، ثم انتقل إلى الخطوة الزمنية التالية.

موضوعات ذات صلة