Langkau tajuk talian

Kaedah Newton-Raphson

Pelinearan dan Hubungan Ulangan Iteratif

Kerja Maya Daya Luaran dan Pemasangan Persamaan Global memberikan persamaan tak linear bagi sesaran nod pada masa \(t_{n+1}\), \(\boldsymbol{u}_{n+1}\), yang diselesaikan dengan kaedah Newton-Raphson. Sesaran nod sehingga masa \(t_n\), \(\boldsymbol{u}_n\), dianggap diketahui, dan kenaikan sesaran \(\Delta\boldsymbol{u}\) diambil sebagai pemboleh ubah tak diketahui yang perlu ditentukan

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

Selepas ini, kebergantungan vektor daya luaran pada sesaran nod 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 penyelesaian semasa \(\Delta\boldsymbol{u}\), takrifkan kekakuan tangen

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

Dengan menggunakan ini, pelinearan persamaan tak linear memberikan

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

Andaikan pembetulan pada lelaran ke-\(i\) ialah \(d\boldsymbol{u}_i\), dan vektor baki pada permulaan lelaran ialah

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

Maka hubungan ulangan iteratif ialah

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

Oleh itu, baki \(\boldsymbol{R}_i\) ialah kuantiti yang sepadan dengan ketakseimbangan daya daripada keadaan keseimbangan.

Pembinaan Matriks Kekakuan Tangen

Kekakuan tangen \(\boldsymbol{K} = \partial\boldsymbol{Q}/\partial\boldsymbol{u}\) dibina dengan membeza separa vektor daya dalaman unsur yang diperoleh dalam Pendiskretan Kerja Maya Daya Dalaman terhadap sesaran nod, mengamirkan integran aras unsur yang terhasil pada setiap domain unsur, dan memasangnya. Dengan melambangkan integran aras unsur sebagai \(\boldsymbol{K}^e_X\) (notasi konfigurasi rujukan, formulasi TL) atau \(\boldsymbol{K}^e_x\) (notasi konfigurasi semasa, formulasi UL), kekakuan tangen unsur ialah

\[ \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. Dalam kedua-dua kes, integran diuraikan kepada jumlah sebutan kekakuan bahan (sebutan sesaran awal) dan sebutan kekakuan geometri (sebutan tegasan awal).

Formulasi Total Lagrange

Dalam formulasi Total Lagrange, hubungan linear diandaikan antara kadar tegasan Piola-Kirchhoff kedua \(\dot{\boldsymbol{S}}\) dan kadar terikan Green-Lagrange \(\dot{\boldsymbol{E}}\), iaitu \(\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}\). Ini sepadan dengan hukum konstitutif bagi bahan elastik linear (bahan St. Venant-Kirchhoff) dan bahan hiperelastik, dan FrontISTR menggunakan formulasi Total Lagrange untuk bahan-bahan ini. Integran kekakuan tangen unsur 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}}) \]

Sebutan pertama di sebelah kanan ialah sebutan kekakuan bahan (sebutan sesaran awal), dan sebutan kedua ialah sebutan kekakuan geometri (sebutan tegasan awal).

Dalam pelaksanaan FrontISTR, integran ini dinilai 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 seperti berikut. \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}\) ialah matriks B yang diperkenalkan dalam Pendiskretan Kerja Maya Daya Dalaman, dan \(\tilde{\boldsymbol{C}}\) ialah perwakilan Voigt bagi tensor konstitutif \(\boldsymbol{\mathsf{C}}\), iaitu matriks kekakuan bahan \(6\times 6\) (Notasi Tensor dan Asas Matematik). \(\boldsymbol{S}_9, \boldsymbol{F}_9\) ialah matriks penyusunan semula berikut yang digunakan untuk menyatakan sebutan kekakuan geometri sebagai hasil darab matriks. Mula-mula, bagi tensor tertib kedua \(3\times 3\) \(\boldsymbol{A}\), takrifkan notasi \([\,\cdot\,]\) yang menyusunnya semula 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 takrif ini, \(\boldsymbol{F}_9\) menyatakan variasi kecerunan ubah bentuk dalam bentuk \([\delta\boldsymbol{F}] = \boldsymbol{F}_9\, \delta\boldsymbol{u}^e\) dan merupakan matriks \(9\times d n_e\). Bagi nod unsur \(\alpha = 1, \ldots, n_e\), blok \(9\times d\) yang sepadan ialah

\[ [\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{ ialah matriks identiti } 3\times 3) \]

dan diberikan oleh \(\boldsymbol{F}_9 = [[\boldsymbol{F}_9]_1, \ldots, [\boldsymbol{F}_9]_{n_e}]\), dengan blok disusun secara mendatar mengikut tertib nod unsur. \(\boldsymbol{S}_9\) dipilih supaya, apabila digabungkan dengan matriks ini, sebutan kekakuan geometri dinyatakan sebagai \(\delta\boldsymbol{u}^{eT}\, \boldsymbol{F}_9^T \boldsymbol{S}_9 \boldsymbol{F}_9\, \dot{\boldsymbol{u}}^e\); ia ialah 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 terhasil.

Formulasi Updated Lagrange

Dalam formulasi Updated Lagrange, hubungan linear diandaikan antara kadar Jaumann bagi tensor tegasan Kirchhoff relatif \(\hat{\boldsymbol{\sigma}}^{\nabla J}\) dan tensor kadar ubah bentuk \(\boldsymbol{D}\), iaitu \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\). Ini ialah bentuk hukum konstitutif hipoelastik yang umum bagi bahan elastik linear, elastoplastik dan rayapan, dan FrontISTR menggunakan formulasi Updated Lagrange untuk bahan-bahan ini. Integran kekakuan tangen unsur yang dinyatakan dalam konfigurasi semasa 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}) \]

dengan \(\boldsymbol{\sigma}^{\nabla T}\) ialah kadar Truesdell, \(\boldsymbol{A}_{(L)}\) ialah bahagian linear terikan Almansi, \(\boldsymbol{F}_t = \partial\boldsymbol{u}/\partial\boldsymbol{x}\) ialah kecerunan sesaran terhadap konfigurasi semasa, dan \(\boldsymbol{L}\) ialah tensor kecerunan halaju. Sebutan pertama di sebelah kanan ialah sebutan kekakuan bahan, dan sebutan kedua ialah sebutan kekakuan geometri.

Dalam pelaksanaan FrontISTR, integran ini dinilai 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}\) ialah matriks B yang dibina dalam konfigurasi semasa (Pendiskretan Kerja Maya Daya Dalaman). \(\boldsymbol{\sigma}_9, \boldsymbol{f}_9\) diperoleh daripada \(\boldsymbol{S}_9, \boldsymbol{F}_9\) yang ditakrifkan untuk formulasi TL dengan menggantikan tegasan PK kedua \(\boldsymbol{S}\) dengan tegasan Cauchy \(\boldsymbol{\sigma}\) dan kecerunan konfigurasi rujukan \(\partial N_\alpha^e/\partial X_i\) dengan kecerunan konfigurasi semasa \(\partial N_\alpha^e/\partial x_i\).

\(\boldsymbol{G}\) ialah matriks pembetulan bergantung tegasan Cauchy yang diperlukan untuk menjadikan hukum konstitutif hipoelastik \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) konsisten dengan rangka kerja kekakuan tangen sebagai hukum konstitutif berasaskan kadar Truesdell. Ia diperoleh dengan menyusun komponen tensor tertib 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 terhasil.

Pemasangan Matriks Kekakuan Global

Kekakuan tangen global \(\boldsymbol{K}\) diperoleh dengan membahagikan setiap kekakuan unsur \(\boldsymbol{K}^e\) kepada blok \(d\times d\) \(\boldsymbol{K}^e_{\alpha\beta}\) bagi setiap pasangan nod dan menggunakan set pemasangan tensor tertib kedua \(\mathcal{E}^2(i_g, i_h)\) yang diperkenalkan dalam Pemasangan Kuantiti Fizik Nod Unsur:

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

Nilai yang terhasil disusun sebagai matriks dengan baris \(i_g\) dan lajur \(i_h\). Dalam pelaksanaan, set \(\mathcal{E}^2\) tidak dibina secara eksplisit; sebaliknya, blok yang sepadan ditambah secara langsung dalam gelung unsur. Matriks itu berbentuk segi empat sama dengan dimensi darjah kebebasan per nod \(\times\) jumlah nod \(n_g\), tetapi kerana komponen selain antara nod yang bersambung melalui unsur ialah \(0\), ia disimpan dalam bentuk matriks jarang.

Matriks kekakuan unsur bagi formulasi TL dan UL mempunyai bentuk yang sama kecuali pertukaran konfigurasi rujukan (koordinat nod dan sumber yang digunakan untuk membina matriks B) serta kewujudan atau ketiadaan matriks \(\boldsymbol{G}\). Oleh itu, FrontISTR melaksanakan kedua-dua formulasi dalam subrutin yang sama.

Algoritma Lelaran

Merumuskan perkara di atas, pada permulaan lelaran tetapkan \(\Delta\boldsymbol{u} = \boldsymbol{0}\) dan kira baki awal \(\boldsymbol{R}_0 = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n)\). Kemudian, pada lelaran ke-\(i\), lakukan prosedur berikut.

  1. Pada sesaran semasa \(\boldsymbol{u}_n + \Delta\boldsymbol{u}\), kira kekakuan tangen \(\boldsymbol{K}_i = \partial\boldsymbol{Q}/\partial\boldsymbol{u}|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}}\) menggunakan prosedur dalam Pembinaan Matriks Kekakuan Tangen.
  2. Untuk mengenakan syarat sempadan geometri, ubah suai matriks kekakuan tangen dan vektor baki bagi darjah kebebasan yang tertakluk pada kekangan sesaran, lalu memperoleh \(\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}\) (Pengendalian Syarat Sempadan Geometri).
  3. Selesaikan persamaan linear \(\tilde{\boldsymbol{K}}_i\, d\boldsymbol{u}_i = \tilde{\boldsymbol{R}}_{i-1}\) untuk mendapatkan pembetulan \(d\boldsymbol{u}_i\). Prosedur ini sering merangkumi sebahagian besar kos pengiraan bagi pengiraan iteratif.
  4. Kemas kini kenaikan sesaran sebagai \(\Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i\), dan seterusnya kira vektor daya dalaman \(\boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\) serta baki \(\boldsymbol{R}_i = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\).
  5. Periksa penumpuan dan tamatkan lelaran jika penumpuan dicapai. Komponen yang sepadan dengan tindak balas terkekang muncul dalam baki \(\boldsymbol{R}_i\) pada darjah kebebasan yang tertakluk pada syarat sempadan geometri, maka penunjuk penumpuan dibina daripada \(\tilde{\boldsymbol{R}}_i\) selepas mengecualikan komponen tersebut. Penunjuk dan ambang penumpuan khusus diterangkan dalam Kriteria Penumpuan. Jika penumpuan tidak dicapai dan had lelaran dicapai, lelaran dianggap gagal.

Apabila lelaran menumpu, tambahkan \(\Delta\boldsymbol{u}\) yang telah menumpu kepada \(\boldsymbol{u}_n\) untuk mendapatkan sesaran terkumpul pada masa \(t_{n+1}\), \(\boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u}\), dan teruskan ke langkah masa berikutnya.

Topik Berkaitan