Bỏ qua

Phương pháp Newton-Raphson

Tuyến tính hóa và quan hệ truy hồi lặp

Công ảo của ngoại lực và lắp ráp phương trình toàn cục cho một phương trình phi tuyến đối với chuyển vị nút tại thời điểm \(t_{n+1}\), \(\boldsymbol{u}_{n+1}\), được giải bằng phương pháp Newton-Raphson. Chuyển vị nút đến thời điểm \(t_n\), \(\boldsymbol{u}_n\), được giả sử đã biết, và gia số chuyển vị \(\Delta\boldsymbol{u}\) được lấy làm biến chưa biết cần xác định

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

Từ đây về sau, bỏ qua sự phụ thuộc của vectơ ngoại lực vào chuyển vị nút, và với \(\boldsymbol{F}(\boldsymbol{u}_{n+1}) = \boldsymbol{F}_{n+1}\),

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

được giải.

Tại nghiệm hiện tại \(\Delta\boldsymbol{u}\), định nghĩa độ cứng tiếp tuyến

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

Dùng đại lượng này, việc tuyến tính hóa phương trình phi tuyến cho

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

Gọi lượng hiệu chỉnh ở lần lặp thứ \(i\) là \(d\boldsymbol{u}_i\), và vectơ phần dư tại đầu lần lặp là

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

Khi đó quan hệ truy hồi lặp là

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

Do đó, phần dư \(\boldsymbol{R}_i\) là đại lượng tương ứng với độ mất cân bằng lực so với trạng thái cân bằng.

Xây dựng ma trận độ cứng tiếp tuyến

Độ cứng tiếp tuyến \(\boldsymbol{K} = \partial\boldsymbol{Q}/\partial\boldsymbol{u}\) được xây dựng bằng cách lấy đạo hàm riêng của vectơ nội lực phần tử thu được trong Rời rạc hóa công ảo của nội lực theo chuyển vị nút, tích phân các biểu thức dưới dấu tích phân ở cấp phần tử trên từng miền phần tử, rồi lắp ráp chúng. Ký hiệu biểu thức dưới dấu tích phân cấp phần tử là \(\boldsymbol{K}^e_X\) (ký hiệu theo cấu hình tham chiếu, công thức TL) hoặc \(\boldsymbol{K}^e_x\) (ký hiệu theo cấu hình hiện tại, công thức UL), độ cứng tiếp tuyến phần tử là

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

Sau đây là dạng cuối cùng của các biểu thức dưới dấu tích phân TL/UL. Trong cả hai trường hợp, chúng được phân tách thành tổng của hạng độ cứng vật liệu (hạng chuyển vị ban đầu) và hạng độ cứng hình học (hạng ứng suất ban đầu).

Công thức Total Lagrange

Trong công thức Total Lagrange, giả sử có quan hệ tuyến tính giữa tốc độ ứng suất Piola-Kirchhoff thứ hai \(\dot{\boldsymbol{S}}\) và tốc độ biến dạng Green-Lagrange \(\dot{\boldsymbol{E}}\), cụ thể \(\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}\). Quan hệ này tương ứng với các định luật cấu thành cho vật liệu đàn hồi tuyến tính (vật liệu St. Venant-Kirchhoff) và vật liệu siêu đàn hồi, và FrontISTR sử dụng công thức Total Lagrange cho các vật liệu này. Khi đó biểu thức dưới dấu tích phân của độ cứng tiếp tuyến phần tử được viết ở dạng tensor là

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

Hạng thứ nhất ở vế phải là hạng độ cứng vật liệu (hạng chuyển vị ban đầu), và hạng thứ hai là hạng độ cứng hình học (hạng ứng suất ban đầu).

Trong cài đặt FrontISTR, biểu thức dưới dấu tích phân này được tính ở dạng ma trận bằng ký hiệu 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 \]

Mỗi ma trận được định nghĩa như sau. \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}\) là các ma trận B được giới thiệu trong Rời rạc hóa công ảo của nội lực, còn \(\tilde{\boldsymbol{C}}\) là biểu diễn Voigt của tensor cấu thành \(\boldsymbol{\mathsf{C}}\), tức ma trận độ cứng vật liệu \(6\times 6\) (Ký hiệu tensor và cơ sở toán học). \(\boldsymbol{S}_9, \boldsymbol{F}_9\) là các ma trận sắp xếp lại sau đây dùng để biểu diễn hạng độ cứng hình học dưới dạng tích ma trận. Trước hết, đối với tensor bậc hai \(3\times 3\) \(\boldsymbol{A}\), định nghĩa ký hiệu \([\,\cdot\,]\) để sắp xếp nó thành vectơ 9 thành phần như sau

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

Với định nghĩa này, \(\boldsymbol{F}_9\) biểu diễn biến phân của gradient biến dạng ở dạng \([\delta\boldsymbol{F}] = \boldsymbol{F}_9\, \delta\boldsymbol{u}^e\) và là ma trận \(9\times d n_e\). Với nút phần tử \(\alpha = 1, \ldots, n_e\), khối \(9\times d\) tương ứng là

\[ [\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{ là ma trận đơn vị } 3\times 3) \]

và được cho bởi \(\boldsymbol{F}_9 = [[\boldsymbol{F}_9]_1, \ldots, [\boldsymbol{F}_9]_{n_e}]\), với các khối được bố trí theo phương ngang theo thứ tự nút phần tử. \(\boldsymbol{S}_9\) được chọn sao cho khi kết hợp với ma trận này, hạng độ cứng hình học được biểu diễn là \(\delta\boldsymbol{u}^{eT}\, \boldsymbol{F}_9^T \boldsymbol{S}_9 \boldsymbol{F}_9\, \dot{\boldsymbol{u}}^e\); đó là ma trận \(9\times 9\) sau

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

Đây là ma trận thu được.

Công thức Updated Lagrange

Trong công thức Updated Lagrange, giả sử có quan hệ tuyến tính giữa tốc độ Jaumann của tensor ứng suất Kirchhoff tương đối \(\hat{\boldsymbol{\sigma}}^{\nabla J}\) và tensor tốc độ biến dạng \(\boldsymbol{D}\), cụ thể \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\). Đây là dạng của định luật cấu thành hypoelastic dùng chung cho vật liệu đàn hồi tuyến tính, đàn hồi-dẻo và từ biến, và FrontISTR sử dụng công thức Updated Lagrange cho các vật liệu này. Khi đó biểu thức dưới dấu tích phân của độ cứng tiếp tuyến phần tử trong cấu hình hiện tại được viết ở dạng tensor là

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

trong đó \(\boldsymbol{\sigma}^{\nabla T}\) là tốc độ Truesdell, \(\boldsymbol{A}_{(L)}\) là phần tuyến tính của biến dạng Almansi, \(\boldsymbol{F}_t = \partial\boldsymbol{u}/\partial\boldsymbol{x}\) là gradient chuyển vị theo cấu hình hiện tại, và \(\boldsymbol{L}\) là tensor gradient vận tốc. Hạng thứ nhất ở vế phải là hạng độ cứng vật liệu, và hạng thứ hai là hạng độ cứng hình học.

Trong cài đặt FrontISTR, biểu thức dưới dấu tích phân này được tính ở dạng ma trận bằng ký hiệu 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 \]

Ở đây, \(\boldsymbol{b}\) là ma trận B được xây dựng trong cấu hình hiện tại (Rời rạc hóa công ảo của nội lực). \(\boldsymbol{\sigma}_9, \boldsymbol{f}_9\) được thu được từ \(\boldsymbol{S}_9, \boldsymbol{F}_9\) đã định nghĩa cho công thức TL bằng cách thay ứng suất PK thứ hai \(\boldsymbol{S}\) bằng ứng suất Cauchy \(\boldsymbol{\sigma}\) và gradient theo cấu hình tham chiếu \(\partial N_\alpha^e/\partial X_i\) bằng gradient theo cấu hình hiện tại \(\partial N_\alpha^e/\partial x_i\).

\(\boldsymbol{G}\) là ma trận hiệu chỉnh phụ thuộc ứng suất Cauchy, cần thiết để định luật cấu thành hypoelastic \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) nhất quán với khung độ cứng tiếp tuyến như một định luật cấu thành dựa trên tốc độ Truesdell. Ma trận này thu được bằng cách sắp xếp các thành phần tensor bậc bốn \(G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}\) theo dạng Voigt \(6\times 6\) như sau

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

Đây là ma trận thu được.

Lắp ráp ma trận độ cứng toàn cục

Độ cứng tiếp tuyến toàn cục \(\boldsymbol{K}\) thu được bằng cách chia mỗi độ cứng phần tử \(\boldsymbol{K}^e\) thành các khối \(d\times d\) \(\boldsymbol{K}^e_{\alpha\beta}\) cho từng cặp nút và dùng tập lắp ráp tensor bậc hai \(\mathcal{E}^2(i_g, i_h)\) được giới thiệu trong Lắp ráp các đại lượng vật lý tại nút phần tử:

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

Các giá trị thu được được sắp thành ma trận với hàng \(i_g\) và cột \(i_h\). Trong cài đặt, tập \(\mathcal{E}^2\) không được xây dựng tường minh; thay vào đó, các khối tương ứng được cộng trực tiếp trong vòng lặp phần tử. Ma trận là ma trận vuông có kích thước bằng số bậc tự do trên mỗi nút \(\times\) tổng số nút \(n_g\), nhưng vì các thành phần ngoài những cặp nút được liên kết qua phần tử bằng \(0\), nó được lưu ở dạng ma trận thưa.

Các ma trận độ cứng phần tử của công thức TL và UL có cùng dạng, ngoại trừ việc chuyển đổi cấu hình tham chiếu (tọa độ nút và nguồn dùng để xây dựng ma trận B) và sự có mặt hay vắng mặt của ma trận \(\boldsymbol{G}\). Vì vậy FrontISTR cài đặt cả hai công thức trong một chương trình con dùng chung.

Thuật toán lặp

Tóm tắt các nội dung trên, tại đầu quá trình lặp đặt \(\Delta\boldsymbol{u} = \boldsymbol{0}\) và tính phần dư ban đầu \(\boldsymbol{R}_0 = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n)\). Sau đó, ở lần lặp thứ \(i\), thực hiện quy trình sau.

  1. Tại chuyển vị hiện tại \(\boldsymbol{u}_n + \Delta\boldsymbol{u}\), tính độ cứng tiếp tuyến \(\boldsymbol{K}_i = \partial\boldsymbol{Q}/\partial\boldsymbol{u}|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}}\) theo quy trình trong Xây dựng ma trận độ cứng tiếp tuyến.
  2. Để áp đặt điều kiện biên hình học, sửa đổi ma trận độ cứng tiếp tuyến và vectơ phần dư cho các bậc tự do chịu ràng buộc chuyển vị, thu được \(\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}\) (Xử lý điều kiện biên hình học).
  3. Giải phương trình tuyến tính \(\tilde{\boldsymbol{K}}_i\, d\boldsymbol{u}_i = \tilde{\boldsymbol{R}}_{i-1}\) để thu được lượng hiệu chỉnh \(d\boldsymbol{u}_i\). Quy trình này thường chiếm phần lớn chi phí tính toán của phép lặp.
  4. Cập nhật gia số chuyển vị theo \(\Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i\), rồi tương ứng tính vectơ nội lực \(\boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\) và phần dư \(\boldsymbol{R}_i = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\).
  5. Kiểm tra hội tụ và kết thúc quá trình lặp nếu đã hội tụ. Các thành phần tương ứng với phản lực ràng buộc xuất hiện trong phần dư \(\boldsymbol{R}_i\) tại các bậc tự do chịu điều kiện biên hình học, vì vậy chỉ tiêu hội tụ được xây dựng từ \(\tilde{\boldsymbol{R}}_i\) sau khi loại các thành phần này. Các chỉ tiêu và ngưỡng hội tụ cụ thể được mô tả trong Tiêu chí hội tụ. Nếu không hội tụ và đạt giới hạn số lần lặp thì lần giải được coi là thất bại.

Khi quá trình lặp hội tụ, cộng \(\Delta\boldsymbol{u}\) đã hội tụ vào \(\boldsymbol{u}_n\) để thu được chuyển vị tích lũy tại thời điểm \(t_{n+1}\), \(\boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u}\), rồi chuyển sang bước thời gian tiếp theo.

Chủ đề liên quan

AI-assisted translation May contain errors Official docs Status