Newton-Raphson विधि¶

रैखिकीकरण और पुनरावृत्त पुनरावर्तन¶

बाहरी बलों का आभासी कार्य और वैश्विक समीकरण का संयोजन समय $t_{n+1}$ पर nodal displacement $\boldsymbol{u}_{n+1}$ के लिए एक nonlinear equation देता है, जिसे Newton-Raphson विधि से हल किया जाता है। समय $t_n$ तक का nodal displacement $\boldsymbol{u}_n$ ज्ञात माना जाता है और निर्धारित किए जाने वाले अज्ञात चर के रूप में displacement increment $\Delta\boldsymbol{u}$ लिया जाता है

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

आगे external-force vector की nodal displacement पर निर्भरता उपेक्षित की जाती है, और $\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}$ पर tangent stiffness परिभाषित करें

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

इसका उपयोग करने पर nonlinear equation का linearization देता है

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

$i$-वीं iteration में correction को $d\boldsymbol{u}_i$ और iteration की शुरुआत में residual vector को

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

मानें। तब iterative recurrence है

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

इस प्रकार residual $\boldsymbol{R}_i$ संतुलन से force imbalance के अनुरूप एक quantity है।

टैन्जेंट कठोरता मैट्रिक्स का निर्माण¶

tangent stiffness $\boldsymbol{K} = \partial\boldsymbol{Q}/\partial\boldsymbol{u}$ को आंतरिक बलों के आभासी कार्य का विविक्तीकरण में प्राप्त element internal-force vector को nodal displacement के सापेक्ष आंशिक अवकलित करके, प्राप्त element-level integrands को प्रत्येक element domain पर integrate करके और assemble करके बनाया जाता है। element-level integrand को $\boldsymbol{K}^e_X$ (reference-configuration notation, TL formulation) या $\boldsymbol{K}^e_x$ (current-configuration notation, UL formulation) से दर्शाने पर element tangent stiffness है

\[ \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 integrands के अंतिम रूप दिए गए हैं। दोनों मामलों में इन्हें material stiffness term (initial-displacement term) और geometric stiffness term (initial-stress term) के योग में विभाजित किया जाता है।

Total Lagrange Formulation¶

Total Lagrange formulation में second Piola-Kirchhoff stress की दर $\dot{\boldsymbol{S}}$ और Green-Lagrange strain rate $\dot{\boldsymbol{E}}$ के बीच linear relation, अर्थात $\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}$, माना जाता है। यह linear elastic materials (St. Venant-Kirchhoff materials) और hyperelastic materials के constitutive laws के अनुरूप है, और FrontISTR इन materials के लिए Total Lagrange formulation उपयोग करता है। तब element tangent-stiffness integrand tensor form में इस प्रकार लिखा जाता है

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

दाएँ पक्ष का पहला term material stiffness term (initial-displacement term) और दूसरा term geometric stiffness term (initial-stress term) है।

FrontISTR implementation में इस integrand का मूल्यांकन Voigt notation का उपयोग करके matrix form में किया जाता है:

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

प्रत्येक matrix निम्न है। $\boldsymbol{B}_L, \boldsymbol{B}_{NL}$ आंतरिक बलों के आभासी कार्य का विविक्तीकरण में प्रस्तुत B matrices हैं, और $\tilde{\boldsymbol{C}}$ constitutive tensor $\boldsymbol{\mathsf{C}}$ का Voigt representation है, अर्थात $6\times 6$ material stiffness matrix (Tensor notation और गणितीय आधार)। $\boldsymbol{S}_9, \boldsymbol{F}_9$ geometric stiffness term को matrix product के रूप में व्यक्त करने के लिए निम्न rearrangement matrices हैं। पहले, $3\times 3$ second-order tensor $\boldsymbol{A}$ के लिए notation $[\,\cdot\,]$ परिभाषित करें जो इसे 9-component vector में rearrange करता है

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

इस परिभाषा के साथ $\boldsymbol{F}_9$ deformation gradient की variation को $[\delta\boldsymbol{F}] = \boldsymbol{F}_9\, \delta\boldsymbol{u}^e$ के रूप में व्यक्त करता है और यह $9\times d n_e$ matrix है। element node $\alpha = 1, \ldots, n_e$ के लिए संबंधित $9\times d$ block है

\[ [\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{ identity matrix है}) \]

और $\boldsymbol{F}_9 = [[\boldsymbol{F}_9]_1, \ldots, [\boldsymbol{F}_9]_{n_e}]$ से दिया जाता है, जहाँ blocks element-node order में horizontally arranged हैं। $\boldsymbol{S}_9$ ऐसा चुना जाता है कि इस matrix के साथ geometric stiffness term $\delta\boldsymbol{u}^{eT}\, \boldsymbol{F}_9^T \boldsymbol{S}_9 \boldsymbol{F}_9\, \dot{\boldsymbol{u}}^e$ के रूप में व्यक्त हो; यह निम्न $9\times 9$ matrix है

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

यह परिणामी matrix है।

Updated Lagrange Formulation¶

Updated Lagrange formulation में relative Kirchhoff stress tensor की Jaumann rate $\hat{\boldsymbol{\sigma}}^{\nabla J}$ और rate-of-deformation tensor $\boldsymbol{D}$ के बीच linear relation, अर्थात $\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}$, माना जाता है। यह linear elastic, elastoplastic और creep materials के लिए सामान्य hypoelastic constitutive law का रूप है, और FrontISTR इन materials के लिए Updated Lagrange formulation उपयोग करता है। current configuration में व्यक्त element tangent-stiffness integrand तब tensor form में इस प्रकार लिखा जाता है

\[ \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 rate है, $\boldsymbol{A}_{(L)}$ Almansi strain का linear part है, $\boldsymbol{F}_t = \partial\boldsymbol{u}/\partial\boldsymbol{x}$ current configuration के सापेक्ष displacement gradient है और $\boldsymbol{L}$ velocity-gradient tensor है। दाएँ पक्ष का पहला term material stiffness term और दूसरा geometric stiffness term है।

FrontISTR implementation में इस integrand का मूल्यांकन Voigt notation का उपयोग करके matrix form में किया जाता है:

\[ \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}$ current configuration में निर्मित B matrix है (आंतरिक बलों के आभासी कार्य का विविक्तीकरण)। $\boldsymbol{\sigma}_9, \boldsymbol{f}_9$ TL formulation के लिए परिभाषित $\boldsymbol{S}_9, \boldsymbol{F}_9$ से second PK stress $\boldsymbol{S}$ को Cauchy stress $\boldsymbol{\sigma}$ तथा reference-configuration gradient $\partial N_\alpha^e/\partial X_i$ को current-configuration gradient $\partial N_\alpha^e/\partial x_i$ से बदलकर प्राप्त होते हैं।

$\boldsymbol{G}$ Cauchy-stress-dependent correction matrix है, जो hypoelastic constitutive law $\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}$ को Truesdell-rate-based constitutive law के रूप में tangent-stiffness framework के साथ consistent बनाने के लिए आवश्यक है। इसे fourth-order tensor components $G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}$ को $6\times 6$ Voigt form में व्यवस्थित करके प्राप्त किया जाता है

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

यह परिणामी matrix है।

वैश्विक कठोरता मैट्रिक्स का संयोजन¶

global tangent stiffness $\boldsymbol{K}$ प्रत्येक element stiffness $\boldsymbol{K}^e$ को node pairs के लिए $d\times d$ blocks $\boldsymbol{K}^e_{\alpha\beta}$ में विभाजित करके और तत्व-नोडल भौतिक राशियों का संयोजन में प्रस्तुत second-order-tensor assembly set $\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} \]

प्राप्त मानों को row $i_g$ और column $i_h$ वाली matrix के रूप में व्यवस्थित किया जाता है। implementation में set $\mathcal{E}^2$ स्पष्ट रूप से निर्मित नहीं किया जाता; इसके बजाय संबंधित blocks element loop के भीतर सीधे जोड़े जाते हैं। matrix square है, जिसकी dimension प्रति node degrees of freedom $\times$ total nodes $n_g$ के बराबर है, लेकिन elements के माध्यम से जुड़े nodes के बीच के components के अलावा अन्य सभी $0$ होने के कारण इसे sparse-matrix form में stored किया जाता है।

TL और UL formulations के element stiffness matrices का रूप समान है, केवल reference configuration (nodal coordinates और B matrix बनाने में प्रयुक्त source) का switch तथा $\boldsymbol{G}$ matrix की उपस्थिति/अनुपस्थिति अलग है। इसलिए FrontISTR दोनों formulations को common subroutine में implement करता है।

पुनरावृत्ति एल्गोरिथ्म¶

ऊपर का सारांश देते हुए, iteration की शुरुआत में $\Delta\boldsymbol{u} = \boldsymbol{0}$ सेट करें और initial residual $\boldsymbol{R}_0 = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n)$ गणना करें। फिर $i$-वीं iteration में निम्न procedure करें।

वर्तमान displacement $\boldsymbol{u}_n + \Delta\boldsymbol{u}$ पर टैन्जेंट कठोरता मैट्रिक्स का निर्माण की procedure से tangent stiffness $\boldsymbol{K}_i = \partial\boldsymbol{Q}/\partial\boldsymbol{u}|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}}$ गणना करें।
geometric boundary conditions लागू करने के लिए displacement constraints के अधीन degrees of freedom के लिए tangent stiffness matrix और residual vector को modify करें और $\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}$ प्राप्त करें (ज्यामितीय सीमा शर्तों का प्रबंधन)।
correction $d\boldsymbol{u}_i$ प्राप्त करने के लिए linear equation $\tilde{\boldsymbol{K}}_i\, d\boldsymbol{u}_i = \tilde{\boldsymbol{R}}_{i-1}$ हल करें। यह procedure अक्सर iterative calculation की अधिकांश computational cost के लिए उत्तरदायी होती है।
displacement increment को $\Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i$ के रूप में update करें, और तदनुसार internal-force vector $\boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})$ तथा residual $\boldsymbol{R}_i = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})$ गणना करें।
convergence जाँचें और converge होने पर iteration समाप्त करें। geometric boundary conditions के अधीन degrees of freedom पर constrained reactions के अनुरूप components residual $\boldsymbol{R}_i$ में दिखाई देते हैं, इसलिए इन components को हटाने के बाद $\tilde{\boldsymbol{R}}_i$ से convergence indicator बनाया जाता है। विशिष्ट convergence indicators और thresholds अभिसरण मानदंड में वर्णित हैं। converge न होने और iteration limit पहुँचने पर iteration failed मानी जाती है।

iteration converge होने पर converged $\Delta\boldsymbol{u}$ को $\boldsymbol{u}_n$ में जोड़कर समय $t_{n+1}$ का accumulated displacement $\boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u}$ प्राप्त करें और अगले time step पर जाएँ।