Lewati ke isi

Metode Newton-Raphson

Linearisasi dan Relasi Rekurensi Iteratif

Kerja Virtual Gaya Eksternal dan Perakitan Persamaan Global memberikan persamaan nonlinear untuk perpindahan nodal pada waktu \(t_{n+1}\), \(\boldsymbol{u}_{n+1}\), yang diselesaikan dengan metode Newton-Raphson. Perpindahan nodal hingga waktu \(t_n\), \(\boldsymbol{u}_n\), diasumsikan diketahui, dan inkremen perpindahan \(\Delta\boldsymbol{u}\) diambil sebagai variabel tak diketahui yang harus ditentukan

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

Selanjutnya, ketergantungan vektor gaya eksternal pada perpindahan nodal diabaikan, dan dengan \(\boldsymbol{F}(\boldsymbol{u}_{n+1}) = \boldsymbol{F}_{n+1}\),

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

diselesaikan.

Pada solusi saat ini \(\Delta\boldsymbol{u}\), definisikan kekakuan tangen

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

Dengan menggunakan ini, linearisasi persamaan nonlinear menghasilkan

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

Misalkan koreksi pada iterasi ke-\(i\) adalah \(d\boldsymbol{u}_i\), dan vektor residual pada awal iterasi adalah

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

Maka relasi rekurensi iteratif adalah

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

Dengan demikian, residual \(\boldsymbol{R}_i\) adalah besaran yang berkorespondensi dengan ketidakseimbangan gaya terhadap keadaan setimbang.

Pembentukan Matriks Kekakuan Tangen

Kekakuan tangen \(\boldsymbol{K} = \partial\boldsymbol{Q}/\partial\boldsymbol{u}\) dibentuk dengan mendiferensialkan secara parsial vektor gaya internal elemen yang diperoleh dalam Diskretisasi Kerja Virtual Gaya Internal terhadap perpindahan nodal, mengintegrasikan integran tingkat elemen yang dihasilkan pada setiap domain elemen, lalu merakitnya. Dengan menyatakan integran tingkat elemen sebagai \(\boldsymbol{K}^e_X\) (notasi konfigurasi referensi, formulasi TL) atau \(\boldsymbol{K}^e_x\) (notasi konfigurasi saat ini, formulasi UL), kekakuan tangen elemen adalah

\[ \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}) \]

Berikut diberikan bentuk akhir integran TL/UL. Pada kedua kasus, integran didekomposisi menjadi jumlah suku kekakuan material (suku perpindahan awal) dan suku kekakuan geometris (suku tegangan awal).

Formulasi Total Lagrange

Dalam formulasi Total Lagrange, diasumsikan hubungan linear antara laju tegangan Piola-Kirchhoff kedua \(\dot{\boldsymbol{S}}\) dan laju regangan Green-Lagrange \(\dot{\boldsymbol{E}}\), yaitu \(\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}\). Ini bersesuaian dengan hukum konstitutif untuk material elastik linear (material St. Venant-Kirchhoff) dan material hiperelastik, dan FrontISTR menggunakan formulasi Total Lagrange untuk material tersebut. Integran kekakuan tangen elemen kemudian ditulis dalam bentuk tensor sebagai

\[ \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}}) \]

Suku pertama pada ruas kanan adalah suku kekakuan material (suku perpindahan awal), dan suku kedua adalah suku kekakuan geometris (suku tegangan awal).

Dalam implementasi FrontISTR, integran ini dievaluasi dalam bentuk matriks menggunakan notasi 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 \]

Setiap matriks adalah sebagai berikut. \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}\) adalah matriks B yang diperkenalkan dalam Diskretisasi Kerja Virtual Gaya Internal, dan \(\tilde{\boldsymbol{C}}\) adalah representasi Voigt dari tensor konstitutif \(\boldsymbol{\mathsf{C}}\), yaitu matriks kekakuan material \(6\times 6\) (Notasi Tensor dan Dasar Matematis). \(\boldsymbol{S}_9, \boldsymbol{F}_9\) adalah matriks pengaturan ulang berikut yang digunakan untuk menyatakan suku kekakuan geometris sebagai perkalian matriks. Pertama, untuk tensor orde kedua \(3\times 3\) \(\boldsymbol{A}\), definisikan notasi \([\,\cdot\,]\) yang mengaturnya ulang menjadi vektor 9-komponen sebagai

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

Dengan definisi ini, \(\boldsymbol{F}_9\) menyatakan variasi gradien deformasi dalam bentuk \([\delta\boldsymbol{F}] = \boldsymbol{F}_9\, \delta\boldsymbol{u}^e\) dan merupakan matriks \(9\times d n_e\). Untuk node elemen \(\alpha = 1, \ldots, n_e\), blok \(9\times d\) yang bersesuaian adalah

\[ [\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{ adalah } 3\times 3 \text{ matriks identitas}) \]

dan diberikan oleh \(\boldsymbol{F}_9 = [[\boldsymbol{F}_9]_1, \ldots, [\boldsymbol{F}_9]_{n_e}]\), dengan blok disusun horizontal menurut urutan node elemen. \(\boldsymbol{S}_9\) dipilih sehingga, ketika digabungkan dengan matriks ini, suku kekakuan geometris dinyatakan sebagai \(\delta\boldsymbol{u}^{eT}\, \boldsymbol{F}_9^T \boldsymbol{S}_9 \boldsymbol{F}_9\, \dot{\boldsymbol{u}}^e\); matriksnya adalah matriks \(9\times 9\) berikut

\[ \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} \]

Inilah matriks yang dihasilkan.

Formulasi Updated Lagrange

Dalam formulasi Updated Lagrange, diasumsikan hubungan linear antara laju Jaumann dari tensor tegangan Kirchhoff relatif \(\hat{\boldsymbol{\sigma}}^{\nabla J}\) dan tensor laju deformasi \(\boldsymbol{D}\), yaitu \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\). Ini merupakan bentuk hukum konstitutif hipoelastik yang umum untuk material elastik linear, elastoplastik, dan creep, dan FrontISTR menggunakan formulasi Updated Lagrange untuk material tersebut. Integran kekakuan tangen elemen yang dinyatakan pada konfigurasi saat ini kemudian ditulis dalam bentuk tensor sebagai

\[ \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}) \]

di mana \(\boldsymbol{\sigma}^{\nabla T}\) adalah laju Truesdell, \(\boldsymbol{A}_{(L)}\) adalah bagian linear regangan Almansi, \(\boldsymbol{F}_t = \partial\boldsymbol{u}/\partial\boldsymbol{x}\) adalah gradien perpindahan terhadap konfigurasi saat ini, dan \(\boldsymbol{L}\) adalah tensor gradien kecepatan. Suku pertama pada ruas kanan adalah suku kekakuan material, dan suku kedua adalah suku kekakuan geometris.

Dalam implementasi FrontISTR, integran ini dievaluasi dalam bentuk matriks menggunakan notasi 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 \]

Di sini, \(\boldsymbol{b}\) adalah matriks B yang dibentuk pada konfigurasi saat ini (Diskretisasi Kerja Virtual Gaya Internal). \(\boldsymbol{\sigma}_9, \boldsymbol{f}_9\) diperoleh dari \(\boldsymbol{S}_9, \boldsymbol{F}_9\) yang didefinisikan untuk formulasi TL dengan mengganti tegangan PK kedua \(\boldsymbol{S}\) dengan tegangan Cauchy \(\boldsymbol{\sigma}\) dan gradien konfigurasi referensi \(\partial N_\alpha^e/\partial X_i\) dengan gradien konfigurasi saat ini \(\partial N_\alpha^e/\partial x_i\).

\(\boldsymbol{G}\) adalah matriks koreksi yang bergantung pada tegangan Cauchy dan diperlukan agar hukum konstitutif hipoelastik \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) konsisten dengan kerangka kekakuan tangen sebagai hukum konstitutif berbasis laju Truesdell. Matriks ini diperoleh dengan menyusun komponen tensor orde keempat \(G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}\) dalam bentuk Voigt \(6\times 6\) sebagai

\[ \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} \]

Inilah matriks yang dihasilkan.

Perakitan Matriks Kekakuan Global

Kekakuan tangen global \(\boldsymbol{K}\) diperoleh dengan membagi setiap kekakuan elemen \(\boldsymbol{K}^e\) menjadi blok \(d\times d\) \(\boldsymbol{K}^e_{\alpha\beta}\) untuk setiap pasangan node dan menggunakan himpunan perakitan tensor orde kedua \(\mathcal{E}^2(i_g, i_h)\) yang diperkenalkan dalam Perakitan Besaran Fisik Nodal Elemen:

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

Nilai yang dihasilkan disusun sebagai matriks dengan baris \(i_g\) dan kolom \(i_h\). Dalam implementasi, himpunan \(\mathcal{E}^2\) tidak dibentuk secara eksplisit; sebagai gantinya, blok yang bersesuaian ditambahkan langsung di dalam loop elemen. Matriks ini berbentuk persegi dengan dimensi sama dengan derajat kebebasan per node \(\times\) jumlah total node \(n_g\), tetapi karena komponen selain yang berada di antara node yang terhubung melalui elemen bernilai \(0\), matriks disimpan dalam bentuk matriks sparse.

Matriks kekakuan elemen untuk formulasi TL dan UL memiliki bentuk yang sama kecuali pada pergantian konfigurasi referensi (koordinat nodal dan sumber yang digunakan untuk membentuk matriks B) serta ada atau tidaknya matriks \(\boldsymbol{G}\). Karena itu, FrontISTR mengimplementasikan kedua formulasi dalam subrutin yang sama.

Algoritme Iterasi

Sebagai rangkuman, pada awal iterasi tetapkan \(\Delta\boldsymbol{u} = \boldsymbol{0}\) dan hitung residual awal \(\boldsymbol{R}_0 = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n)\). Kemudian, pada iterasi ke-\(i\), lakukan prosedur berikut.

  1. Pada perpindahan saat ini \(\boldsymbol{u}_n + \Delta\boldsymbol{u}\), hitung kekakuan tangen \(\boldsymbol{K}_i = \partial\boldsymbol{Q}/\partial\boldsymbol{u}|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}}\) menggunakan prosedur dalam Pembentukan Matriks Kekakuan Tangen.
  2. Untuk menerapkan kondisi batas geometris, modifikasi matriks kekakuan tangen dan vektor residual untuk derajat kebebasan yang dikenai kendala perpindahan, sehingga diperoleh \(\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}\) (Penanganan Kondisi Batas Geometris).
  3. Selesaikan persamaan linear \(\tilde{\boldsymbol{K}}_i\, d\boldsymbol{u}_i = \tilde{\boldsymbol{R}}_{i-1}\) untuk memperoleh koreksi \(d\boldsymbol{u}_i\). Prosedur ini sering menyumbang sebagian besar biaya komputasi dari perhitungan iteratif.
  4. Perbarui inkremen perpindahan menjadi \(\Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i\), dan sesuai dengan itu hitung vektor gaya internal \(\boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\) dan residual \(\boldsymbol{R}_i = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\).
  5. Periksa konvergensi dan akhiri iterasi jika konvergensi tercapai. Komponen yang bersesuaian dengan reaksi terkendala muncul dalam residual \(\boldsymbol{R}_i\) pada derajat kebebasan yang dikenai kondisi batas geometris, sehingga indikator konvergensi dibentuk dari \(\tilde{\boldsymbol{R}}_i\) setelah komponen tersebut dikeluarkan. Indikator dan ambang konvergensi spesifik dijelaskan dalam Kriteria Konvergensi. Jika konvergensi tidak tercapai dan batas iterasi tercapai, iterasi dianggap gagal.

Ketika iterasi konvergen, tambahkan \(\Delta\boldsymbol{u}\) yang telah konvergen ke \(\boldsymbol{u}_n\) untuk memperoleh perpindahan kumulatif pada waktu \(t_{n+1}\), \(\boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u}\), lalu lanjutkan ke time step berikutnya.

Topik Terkait