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}\) ला 9 घटकांच्या सदिशात पुनर्मांडणारे \([\,\cdot\,]\) हे चिन्ह

\[ [\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 पद्धतीत सापेक्ष Kirchhoff ताण टेन्सरचा Jaumann दर \(\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\) हे TL पद्धतीसाठी परिभाषित \(\boldsymbol{S}_9, \boldsymbol{F}_9\) मध्ये दुसरा PK ताण \(\boldsymbol{S}\) ऐवजी Cauchy ताण \(\boldsymbol{\sigma}\) आणि संदर्भ विन्यास प्रवणता \(\partial N_\alpha^e/\partial X_i\) ऐवजी सध्याच्या विन्यासाची प्रवणता \(\partial N_\alpha^e/\partial x_i\) वापरून मिळतात.

\(\boldsymbol{G}\) हा हायपोइलॅस्टिक पदार्थ-वर्तन संबंध \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) ला Truesdell-दर-आधारित पदार्थ-वर्तन संबंध म्हणून स्पर्शक कठोरता चौकटीशी सुसंगत करण्यासाठी आवश्यक Cauchy ताणावर अवलंबून असलेला दुरुस्ती मॅट्रिक्स आहे. चौथ्या क्रमाच्या टेन्सरचे घटक \(G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}\) यांना \(6\times 6\) Voigt रूपात मांडल्यास

\[ \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}\) मिळवले जाते आणि पुढील काल-पायरीकडे जाता येते.