შიგთავსზე გადასვლა

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}}\) — კონსტიტუციური ტენზორის \(\boldsymbol{\mathsf{C}}\) Voigt-ის წარმოდგენა, ანუ \(6\times 6\) მასალის სიხისტის მატრიცა (ტენზორული აღნიშვნა და მათემატიკური საფუძვლები). \(\boldsymbol{S}_9, \boldsymbol{F}_9\) არის შემდეგი გადალაგების მატრიცები, რომლებიც გეომეტრიული სიხისტის წევრს მატრიცული ნამრავლით გამოსახავს. ჯერ \(3\times 3\) მეორე რანგის ტენზორისთვის \(\boldsymbol{A}\) განვსაზღვროთ აღნიშვნა \([\,\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-ის ფორმულირებაში დაშვებულია წრფივი კავშირი ფარდობითი Kirchhoff-ის ძაბვის ტენზორის Jaumann-ის სიჩქარეს \(\hat{\boldsymbol{\sigma}}^{\nabla J}\) და დეფორმაციის სიჩქარის ტენზორს \(\boldsymbol{D}\) შორის: \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\). ეს არის წრფივი ელასტიკური, ელასტოპლასტიკური და ცოცვადობის მასალებისთვის საერთო ჰიპოელასტიკური კონსტიტუციური კანონის ფორმა, და 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\) მიიღება TL ფორმულირებისთვის განსაზღვრული \(\boldsymbol{S}_9, \boldsymbol{F}_9\)-იდან მეორე PK ძაბვის \(\boldsymbol{S}\) Cauchy-ის ძაბვით \(\boldsymbol{\sigma}\) და საცნობარო კონფიგურაციის გრადიენტის \(\partial N_\alpha^e/\partial X_i\) მიმდინარე კონფიგურაციის გრადიენტით \(\partial N_\alpha^e/\partial x_i\) ჩანაცვლებით.

\(\boldsymbol{G}\) არის Cauchy-ის ძაბვაზე დამოკიდებული კორექციის მატრიცა, რომელიც საჭიროა ჰიპოელასტიკური კონსტიტუციური კანონის \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) Truesdell-ის სიჩქარეზე დაფუძნებული მხები სიხისტის ჩარჩოსთან შესათავსებლად. იგი მიიღება მეოთხე რანგის ტენზორის კომპონენტების \(G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}\) \(6\times 6\) Voigt ფორმაში დალაგებით:

\[ \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}\), რის შემდეგაც ანალიზი შემდეგ დროით ბიჯზე გადადის.

დაკავშირებული თემები