Ana içeriğe geç

Newton-Raphson Yöntemi

Doğrusallaştırma ve Yineleme Bağıntısı

Dış Kuvvetlerin Sanal İşi ve Global Denklemlerin Montajı bölümünde elde edilen, \(t_{n+1}\) zamanındaki düğüm yer değiştirmesi \(\boldsymbol{u}_{n+1}\) için doğrusal olmayan denklem Newton-Raphson yöntemiyle çözülür. \(t_n\) zamanına kadar düğüm yer değiştirmesi \(\boldsymbol{u}_n\) biliniyor kabul edilir ve yer değiştirme artımı \(\Delta\boldsymbol{u}\) bilinmeyen olarak alınarak

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

elde edilir. Bundan sonra dış kuvvet vektörünün düğüm yer değiştirmesine bağımlılığı ihmal edilir ve \(\boldsymbol{F}(\boldsymbol{u}_{n+1}) = \boldsymbol{F}_{n+1}\) kabul edilerek

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

çözülür.

Mevcut çözüm \(\Delta\boldsymbol{u}\) için teğet rijitlik

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

kullanılarak doğrusal olmayan denklem doğrusallaştırıldığında

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

elde edilir. \(i\). yinelemedeki düzeltme \(d\boldsymbol{u}_i\), yineleme başlangıcındaki artık vektörü ise

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

olarak yazılırsa yineleme bağıntısı

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

biçimini alır. Artık \(\boldsymbol{R}_i\), dengeden sapmaya karşılık gelen kuvvet dengesizliği büyüklüğüdür.

Teğet Rijitlik Matrisinin Oluşturulması

Teğet rijitlik \(\boldsymbol{K} = \partial\boldsymbol{Q}/\partial\boldsymbol{u}\), İç Kuvvetlerin Sanal İşinin Ayrıklaştırılması bölümünde elde edilen eleman iç kuvvet vektörünün düğüm yer değiştirmesine göre kısmi türevi alınarak, ortaya çıkan eleman düzeyi integrandların her eleman bölgesi üzerinde integrasyonu ve montajıyla oluşturulur. Eleman düzeyi integrand \(\boldsymbol{K}^e_X\) (referans konfigürasyon gösterimi, TL yöntemi) veya \(\boldsymbol{K}^e_x\) (mevcut konfigürasyon gösterimi, UL yöntemi) ile gösterildiğinde eleman teğet rijitliği

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

ile verilir. Aşağıda TL/UL integrandlarının son biçimleri gösterilmiştir. Her ikisi de malzeme rijitliği terimi (başlangıç yer değiştirmesi terimi) ile geometrik rijitlik teriminin (başlangıç gerilmesi terimi) toplamına ayrılır.

Total Lagrange Formülasyonu

Total Lagrange formülasyonunda ikinci Piola-Kirchhoff gerilme hızı \(\dot{\boldsymbol{S}}\) ile Green-Lagrange gerinim hızı \(\dot{\boldsymbol{E}}\) arasında \(\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}\) doğrusal ilişkisi varsayılır. Bu, doğrusal elastik (St.Venant-Kirchhoff) ve hiperelastik malzemelerin kurucu bağıntılarına karşılık gelir; FrontISTR Total Lagrange formülasyonunu bu malzemeler için kullanır. Eleman teğet rijitliği integrandı tensör biçiminde

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

olarak yazılır. Sağ taraftaki ilk terim malzeme rijitliği terimi (başlangıç yer değiştirmesi terimi), ikinci terim geometrik rijitlik terimidir (başlangıç gerilmesi terimi).

FrontISTR uygulamasında bu integrand Voigt gösterimi kullanılarak matris biçiminde

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

hesaplanır. Matrisler şöyledir. \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}\), İç Kuvvetlerin Sanal İşinin Ayrıklaştırılması bölümünde tanıtılan B matrisleridir; \(\tilde{\boldsymbol{C}}\), kurucu tensör \(\boldsymbol{\mathsf{C}}\)'nin Voigt gösterimindeki \(6\times 6\) malzeme rijitlik matrisidir (Tensör Gösterimi ve Matematiksel Temeller). \(\boldsymbol{S}_9, \boldsymbol{F}_9\), geometrik rijitlik terimini matris çarpımı olarak ifade etmek için kullanılan yeniden düzenleme matrisleridir. Önce \(3\times 3\) ikinci mertebe tensör \(\boldsymbol{A}\)'yı 9 bileşenli vektöre düzenleyen \([\,\cdot\,]\) gösterimi

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

olarak tanımlandığında, \(\boldsymbol{F}_9\) deformasyon gradyanının varyasyonunu \([\delta\boldsymbol{F}] = \boldsymbol{F}_9\, \delta\boldsymbol{u}^e\) biçiminde ifade eden \(9\times d n_e\) matristir. Eleman düğümü \(\alpha = 1, \ldots, n_e\) için karşılık gelen \(9\times d\) blok

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

olup blokların eleman düğümü sırasıyla yatay dizilmesiyle \(\boldsymbol{F}_9 = [[\boldsymbol{F}_9]_1, \ldots, [\boldsymbol{F}_9]_{n_e}]\) elde edilir. \(\boldsymbol{S}_9\), bununla birlikte geometrik rijitlik teriminin \(\delta\boldsymbol{u}^{eT}\, \boldsymbol{F}_9^T \boldsymbol{S}_9 \boldsymbol{F}_9\, \dot{\boldsymbol{u}}^e\) olarak ifade edilmesini sağlayacak şekilde seçilen aşağıdaki \(9\times 9\) matristir

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

Bu matris elde edilir.

Updated Lagrange Formülasyonu

Updated Lagrange formülasyonunda göreli Kirchhoff gerilme tensörünün Jaumann hızı \(\hat{\boldsymbol{\sigma}}^{\nabla J}\) ile deformasyon hızı tensörü \(\boldsymbol{D}\) arasında \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) doğrusal ilişkisi varsayılır. Bu, doğrusal elastik, elastoplastik ve sünme malzemelerinde ortak hipoelastik kurucu bağıntı biçimidir; FrontISTR Updated Lagrange formülasyonunu bu malzemeler için kullanır. Mevcut konfigürasyondaki eleman teğet rijitliği integrandı tensör biçiminde

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

olarak yazılır (\(\boldsymbol{\sigma}^{\nabla T}\) Truesdell hızı, \(\boldsymbol{A}_{(L)}\) Almansi geriniminin doğrusal kısmı, \(\boldsymbol{F}_t = \partial\boldsymbol{u}/\partial\boldsymbol{x}\) mevcut konfigürasyona göre yer değiştirme gradyanı ve \(\boldsymbol{L}\) hız gradyanı tensörüdür). Sağ taraftaki ilk terim malzeme rijitliği terimi, ikinci terim geometrik rijitlik terimidir.

FrontISTR uygulamasında bu integrand Voigt gösterimi kullanılarak matris biçiminde

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

hesaplanır. \(\boldsymbol{b}\) mevcut konfigürasyonda oluşturulan B matrisidir (İç Kuvvetlerin Sanal İşinin Ayrıklaştırılması). \(\boldsymbol{\sigma}_9, \boldsymbol{f}_9\), TL formülasyonunda tanımlanan \(\boldsymbol{S}_9, \boldsymbol{F}_9\) içinde ikinci PK gerilmesi \(\boldsymbol{S}\) yerine Cauchy gerilmesi \(\boldsymbol{\sigma}\) ve referans-konfigürasyon gradyanı \(\partial N_\alpha^e/\partial X_i\) yerine mevcut-konfigürasyon gradyanı \(\partial N_\alpha^e/\partial x_i\) konularak elde edilir.

\(\boldsymbol{G}\), hipoelastik kurucu bağıntı \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\)'yi Truesdell hızı tabanlı bir kurucu bağıntı olarak teğet rijitlik çerçevesiyle tutarlı kılmak için gereken Cauchy-gerilmesine bağlı düzeltme matrisidir; dördüncü mertebe tensör bileşenleri \(G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}\) \(6\times 6\) Voigt biçiminde

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

Bu matris elde edilir.

Global Rijitlik Matrisinin Montajı

Global teğet rijitlik \(\boldsymbol{K}\), her eleman rijitliği \(\boldsymbol{K}^e\)'nin düğüm çiftleri için \(d\times d\) bloklar \(\boldsymbol{K}^e_{\alpha\beta}\) halinde ayrılması ve Eleman Düğüm Büyüklüklerinin Montajı bölümünde tanıtılan ikinci mertebe tensör montaj kümesi \(\mathcal{E}^2(i_g, i_h)\) kullanılarak

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

elde edilen değerlerin satır \(i_g\), sütun \(i_h\) olacak şekilde düzenlenmesiyle oluşur. Uygulamada \(\mathcal{E}^2\) kümesi açıkça oluşturulmaz; karşılık gelen bloklar eleman döngüsünde doğrudan eklenir. Matris, düğüm başına serbestlik derecesi \(\times\) toplam düğüm sayısı \(n_g\) boyutunda kare matristir; ancak elemanlar üzerinden bağlı olmayan düğümler arasındaki bileşenler \(0\) olduğundan seyrek matris biçiminde saklanır.

TL ve UL formülasyonlarının eleman rijitlikleri, referans konfigürasyonun değiştirilmesi (düğüm koordinatları ve B matrisinin oluşturulduğu kaynak) ve \(\boldsymbol{G}\) matrisinin bulunup bulunmaması dışında aynı biçimdedir. Bu nedenle FrontISTR her ikisini ortak bir alt yordamda uygular.

Yineleme Algoritması

Yukarıdakiler özetlendiğinde, yineleme başlangıcında \(\Delta\boldsymbol{u} = \boldsymbol{0}\) alınır ve başlangıç artığı \(\boldsymbol{R}_0 = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n)\) hesaplanır. Ardından \(i\). yinelemede aşağıdaki adımlar uygulanır.

  1. Mevcut yer değiştirme \(\boldsymbol{u}_n + \Delta\boldsymbol{u}\) için teğet rijitlik \(\boldsymbol{K}_i = \partial\boldsymbol{Q}/\partial\boldsymbol{u}|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}}\), Teğet Rijitlik Matrisinin Oluşturulması prosedürüyle hesaplanır.
  2. Geometrik sınır koşullarını uygulamak için, yer değiştirme kısıtı bulunan serbestlik derecelerinde teğet rijitlik matrisi ve artık vektörü değiştirilerek \(\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}\) elde edilir (Geometrik Sınır Koşullarının İşlenmesi).
  3. \(\tilde{\boldsymbol{K}}_i\, d\boldsymbol{u}_i = \tilde{\boldsymbol{R}}_{i-1}\) doğrusal denklemi çözülerek düzeltme \(d\boldsymbol{u}_i\) elde edilir. Bu işlem çoğu zaman yinelemeli hesabın hesaplama maliyetinin büyük bölümünü oluşturur.
  4. Yer değiştirme artımı \(\Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i\) olarak güncellenir; buna göre iç kuvvet vektörü \(\boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\) ve artık \(\boldsymbol{R}_i = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\) hesaplanır.
  5. Yakınsama kontrol edilir ve sağlanmışsa yineleme sonlandırılır. Geometrik sınır koşulu uygulanan serbestlik derecelerinde artık \(\boldsymbol{R}_i\) içinde kısıt tepki kuvvetine karşılık gelen bileşenler bulunduğundan, yakınsama göstergesi bunlar çıkarıldıktan sonra \(\tilde{\boldsymbol{R}}_i\) bileşenlerinden oluşturulur. Ayrıntılı göstergeler ve eşikler Yakınsama Ölçütleri bölümünde açıklanır. Yakınsama sağlanmadan yineleme üst sınırına ulaşılırsa yineleme başarısız kabul edilir.

Yineleme yakınsadığında, yakınsamış \(\Delta\boldsymbol{u}\) değeri \(\boldsymbol{u}_n\)'ye eklenerek \(t_{n+1}\) zamanındaki birikimli yer değiştirme \(\boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u}\) elde edilir ve sonraki zaman adımına geçilir.

İlgili Konular