วิธี 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 สมมติความสัมพันธ์เชิงเส้นระหว่างอัตรา Jaumann ของเทนเซอร์ความเค้น Kirchhoff สัมพัทธ์ \(\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\) ได้จาก \(\boldsymbol{S}_9, \boldsymbol{F}_9\) ที่นิยามสำหรับสูตร TL โดยแทนความเค้น 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}\) ให้อยู่ในรูป Voigt \(6\times 6\) ดังนี้

\[ \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\) ให้ดำเนินการดังต่อไปนี้

ที่การกระจัดปัจจุบัน \(\boldsymbol{u}_n + \Delta\boldsymbol{u}\) คำนวณสติฟเนสแทนเจนต์ \(\boldsymbol{K}_i = \partial\boldsymbol{Q}/\partial\boldsymbol{u}|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}}\) ตามขั้นตอนใน การสร้างเมทริกซ์สติฟเนสแทนเจนต์
เพื่อบังคับเงื่อนไขขอบเขตเชิงเรขาคณิต ให้แก้ไขเมทริกซ์สติฟเนสแทนเจนต์และเวกเตอร์เศษเหลือสำหรับองศาอิสระที่ถูกจำกัดการกระจัด เพื่อให้ได้ \(\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}\) (การจัดการเงื่อนไขขอบเขตเชิงเรขาคณิต)
แก้สมการเชิงเส้น \(\tilde{\boldsymbol{K}}_i\, d\boldsymbol{u}_i = \tilde{\boldsymbol{R}}_{i-1}\) เพื่อหาค่าปรับแก้ \(d\boldsymbol{u}_i\) ขั้นตอนนี้มักเป็นส่วนที่ใช้ต้นทุนการคำนวณส่วนใหญ่ของการวนซ้ำ
อัปเดตส่วนเพิ่มการกระจัดเป็น \(\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})\) ตามนั้น
ตรวจสอบการลู่เข้าและยุติการวนซ้ำเมื่อได้การลู่เข้า ที่องศาอิสระซึ่งอยู่ภายใต้เงื่อนไขขอบเขตเชิงเรขาคณิต เศษเหลือ \(\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}\) แล้วดำเนินไปยังขั้นเวลาถัดไป

หัวข้อที่เกี่ยวข้อง¶

งานเสมือนของแรงภายนอกและการประกอบสมการส่วนกลาง — จุดเริ่มต้นของสมการไม่เชิงเส้นที่ต้องแก้
การทำให้เป็นดิสครีตของงานเสมือนของแรงภายใน — การสร้าง \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}, \boldsymbol{b}\)
การจัดการเงื่อนไขขอบเขตเชิงเรขาคณิต — การแก้ไขเมทริกซ์สติฟเนสแทนเจนต์และเวกเตอร์เศษเหลือเพื่อบังคับข้อจำกัดการกระจัด
เกณฑ์การลู่เข้า — เกณฑ์หยุดตามนอร์มของเศษเหลือ
สัญกรณ์เทนเซอร์และพื้นฐานทางคณิตศาสตร์ — การแทนเมทริกซ์วัสดุ \(\tilde{\boldsymbol{C}}\) แบบ Voigt
การวนซ้ำไม่เชิงเส้นและการอินทิเกรตเวลา (ฟังก์ชัน) — วิธีใช้และการเลือกในส่วนอ้างอิงฟังก์ชัน