非线性本构关系下圆板的大挠度问题

Large Deflection of Circular Plate under Nonlinear Constitutive Relations

摘要: 本文基于经典板理论和本构非线性关系，研究了周边固支和周边简支圆板在横向均布载荷的作用下的非线性静态响应。假设了材料弹性模量是应变的线性函数，利用能量法推导了本构关系非线性的情况下，圆板弯曲的基本方程和边界条件，对基本方程和边界条件进行了无量纲化处理，并采用打靶方法求解了周边固支和周边简支条件下，圆板轴对称弯曲的数值解。分析了本构非线性参数对圆板轴对称弯曲变形和中性层位置的影响。

Abstract: In this paper, based on the classical plate theory and nonlinear constitutive relation, the nonlinear static response of circular plates with fixed and simply supported edges peripheries under laterally uniform load are studied. It is assumed that the elastic modulus of material is a linear function of strain. The basic equations and boundary conditions of circular plate bending under the condition of non-linear constitutive relation are derived by energy method. The equations and boundary conditions are nondimensionalized, and the numerical solutions of circular plate axisymmetric bending under the conditions of peripheral fixed support and simple support are solved by shooting method. The influence of constitutive nonlinear parameters on bending deformation and neutral layer is analyzed.

1. 引言

由于各项同性及复合材料的力学行为研究在解决工程问题上起到了指导性的作用，国内外科研人员对此进行了大量的研究，并取得了丰硕的成果。但在过去的研究中大多采用经典力学的研究方法，而且多数材料的弹性模量在构件变形过程中变化很小，研究人员对其进行了线性化处理，即在弹性范围内应力应变的本构关系是线性的，而对于本构的非线性关系方向的研究相对较少。实际上材料的本构关系并非完全可近似成线性的。在有些本构非线性分析中采用了双模量(受拉和受压时分段线性)，这一研究理论由苏联学者Ambartsumyan [1] 提出。张鹏、范存新 [2] 研究了在外载荷作用下，四边固支双模量矩形板的弯曲问题，结果表明当材料拉压模量相差较大时，用双模量计算精确度更高，双模量的影响不容忽视。吴晓等 [3] [4] 基于静力平衡方程，确定了拉压不同弹性模量圆板的中性层位置，进一步得出了中心挠度与均布荷载的关系式。并将此结果对比了有限元方法结果。分析表明，拉压弹性模量差别较大的材料在问题求解中，应采用不同拉压弹性模量的弹性理论来求解更精确。近年来，国内外众多学者采用双模量理论对弹性构件的力学行为进行了研究 [5] - [11]，但这些研究都未给出精确的本构非线性关系。而在变形过程中对于非线性本构关系的结构来说，其刚度也随之增加。本文假设了材料的弹性模量是应变的线性函数，并对圆板的轴对称问题进行了分析。

2. 基本方程

本文研究一半径为b，厚度为h，有非线性应力–应变本构关系的圆板。圆板所受横向均布载荷为q。圆柱坐标系(r, φ, z)建立在几何中面上，坐标中心在圆板圆心。根据von Karman圆板大挠度理论，圆板轴对称变形的几何方程为：

${\epsilon }_r=\frac{\text{d}u}{\text{d}r}+\frac{1}{2}{\left(\frac{\text{d}w}{\text{d}r}\right)}^2-z\frac{{\text{d}}^{2}w}{\text{d}{r}^2}$ (1)

${\epsilon }_{\theta }=\frac{u}{r}-z\frac{1}{r}\frac{\text{d}w}{\text{d}r}$ (2)

其中u和w分别为几何中面上点的面内和横向位移。z是厚度方向的坐标 $\left(-h/2\le z\le h/2\right)$ 且规定向下为正。

我们假设弹性模量是应变的线性函数(线性弹性模量理论)即：

$E_r=E\left({\epsilon }_r\right)=E_0+E_1{\epsilon }_r=E_0\left(1+E_P{\epsilon }_r\right)$ (3)

$E_{\theta }=E\left({\epsilon }_{\theta }\right)=E_0+E_1{\epsilon }_{\theta }=E_0\left(1+E_P{\epsilon }_{\theta }\right)$ (4)

其中， $E_0$ 表示材料初始变形时的弹性模量， $E_1$ 表示材料弹性模量随应变的线性变化参数。令 $E_P=\frac{E_1}{E_0}$，该参数体现了本构关系的非线性特性。

由Hooke定律表达式，

(5)

${\epsilon }_{\theta }=\frac{1}{E_{\theta }}\left({\sigma }_{\theta }-\mu {\sigma }_r\right)$ (6)

其中 $E_r$， $E_{\theta }$， $\mu$ 分别为径向和环向的弹性模量及泊松比，其应力场为：

${\sigma }_r=\frac{1}{1-{\mu }^2}\left({\epsilon }_r{E}_r+\mu {\epsilon }_{\theta }{E}_{\theta }\right)$ (7)

${\sigma }_{\theta }=\frac{1}{1-{\mu }^2}\left({\epsilon }_{\theta }{E}_{\theta }+\mu {\epsilon }_r{E}_r\right)$ (8)

内力表达式为：

$\left(N_r,M_r\right)={\displaystyle {\int }_{-h/2}^{h/2}{\sigma }_r\left(1,z\right)\text{d}z}$ (9)

$\left(N_{\theta },M_{\theta }\right)={\displaystyle {\int }_{-h/2}^{h/2}{\sigma }_{\theta }\left(1,z\right)\text{d}z}$ (10)

其中， $N_r$， $N_{\theta }$ 表示面内薄膜力， $M_r$， $M_{\theta }$ 表示弯矩。

圆板变形时的总势能为 $\Pi =U+V$，

其中应变能为：

$U={\displaystyle {\int }_V\left({\displaystyle {\int }_0^{{\epsilon }_r}{\sigma }_r\text{d}\epsilon }+{\displaystyle {\int }_0^{{\epsilon }_{\theta }}{\sigma }_{\theta }\text{d}\epsilon }\right)\text{d}V}$ (11)

外力势能为：

$V=-{\displaystyle {\iint }_{A}qw\text{d}A}+{\displaystyle {\oint }_{\Gamma }pu\text{d}s}$ (12)

根据能量变分原理 $\delta \Pi =\delta U+\delta V=0$ 可得：

$\begin{array}{l}\frac{1}{r}\frac{\text{d}}{\text{d}r}\left(r{N}_r\right)-\frac{N_{\theta }}{r}=0\\ \frac{1}{r}\frac{\text{d}}{\text{d}r}\left(r{N}_r\frac{\text{d}w}{\text{d}r}\right)-\frac{1}{r}\frac{\text{d}{M}_{\theta }}{\text{d}r}+\frac{1}{r}\frac{{\text{d}}^2}{\text{d}{r}^2}\left(r{M}_r\right)+q\end{array}$ (13)

边界条件：

$u=0$ 或 $N_r+p=0$ (14)

$N_r+p=0$ 或 $r{N}_r\frac{\text{d}w}{\text{d}r}-M_{\theta }+r\frac{\text{d}{M}_r}{\text{d}r}=0$ (15)

$\frac{\text{d}w}{\text{d}r}=0$ 或 $M_r=0$ (16)

3. 对基本方程的无量纲化

为了计算方便，而且使求解的结果更具有普遍性，需要对基本方程和边界条件进行无量纲处理，引入下列无量纲参数：

$X=\frac{r}{b},W=\frac{w}{h},U=\frac{bu}{h^2},\delta =\frac{h}{b},Q=\frac{12\left(1-{\mu }^2\right)q}{E_0{\delta }^4}$

对以上基本方程和边界条件进行无量纲化处理，可得无量纲化以后的基本方程和边界条件，不再赘述。

我们采用打靶法数值对问题的无量纲方程进行了求解。令：

$\left\{Y\right\}={\left\{y_1,y_2,y_3,y_4,y_5,y_6\right\}}^{\text{T}}={\left\{U,\frac{\text{d}U}{\text{d}X},W,\frac{\text{d}W}{\text{d}X},\frac{{\text{d}}^{2}W}{\text{d}{X}^2},\frac{{\text{d}}^{3}W}{\text{d}{X}^3}\right\}}^{\text{T}}$

我无量纲化后的基本方程就可化作以下一阶非线性常微分方程：

$\frac{\text{d}Y}{\text{d}X}={\left\{y_2,{\phi }_1,y_4,y_5,y_6,{\phi }_2\right\}}^{\text{T}}$

其中： ${\phi }_1=\frac{s_1}{s_2},{\phi }_2=\frac{s_3+s_4}{s_5}$。

这里由无量纲基本方程可以解得：

$\begin{array}{c}{s}_1=-\frac{12}{X}{y}_2+\frac{6\mu }{X}{y}_4^2-\frac{6}{X}{y}_4^2+\frac{12E_p{\delta }^2\mu }{X}{y}_2^2-\frac{12E_p{\delta }^2}{X}{y}_2^2+\frac{12E_p{\delta }^2\mu }{X}{y}_2{y}_4^2\\ -\frac{12E_p{\delta }^2}{X}{y}_2{y}_4^2+\frac{3E_p{\delta }^2\mu }{X}{y}_4^4-\frac{3E_p{\delta }^2}{X}{y}_4^4+\frac{12}{X^2}{y}_1+\frac{12E_p{\delta }^2\mu }{X^3}{y}_1^2\\ +\frac{12E_p{\delta }^2}{X^3}{y}_1^2-24E_p{\delta }^2{y}_2{y}_4{y}_5-12y_4{y}_5-12E_p{\delta }^2{y}_4^3{y}_5-\frac{24E_p{\delta }^2\mu }{X^2}{y}_1{y}_2\\ +\frac{E_p{\delta }^2\mu }{X^3}{y}_4^2+\frac{E_p{\delta }^2}{X^3}{y}_4^2+\frac{E_p{\delta }^2\mu }{X}{y}_5^2-\frac{E_p{\delta }^2}{X}{y}_5^2-2E_p{\delta }^2{y}_5{y}_6-\frac{2E_p{\delta }^2\mu }{X^2}{y}_4{y}_5\end{array}$

$s_2=12\left(1+2E_p{\delta }^2{y}_2+E_p{\delta }^2{y}_4^2\right)$

$\begin{array}{l}{s}_3=[\frac{24E_p{\delta }^2\mu }{X}{y}_2{y}_4{y}_5-\frac{12}{X}{y}_4{y}_5-\frac{6\mu }{X^2}{y}_4^2-2E_p{\delta }^2{y}_6^2-24E_p{\delta }^2{\phi }_1^2-\frac{24}{X^3}{y}_1+\frac{6}{X^2}{y}_4^2\\ -\frac{12}{X}{\phi }_1-12y_4{y}_6+\frac{24}{X^2}{y}_2-12y_5^2-\frac{3E_p{\delta }^2\mu }{X^2}{y}_4^4-48E_p{\delta }^2{y}_4{y}_5{\phi }_1-\frac{12E_p{\delta }^2}{X}{y}_4^2{\phi }_1\\ -\frac{36E_p{\delta }^2\mu }{X^2}{y}_2^2-24E_p{\delta }^2{y}_2{y}_4{y}_6+\frac{12E_p{\delta }^2}{X^2}{y}_2{y}_4^2-\frac{24E_p{\delta }^2}{X}{y}_2{\phi }_1+\frac{24E_p{\delta }^2}{X^3}{y}_1{y}_2\\ -\frac{36E_p{\delta }^2\mu }{X^4}{y}_1^2-\frac{3E_p{\delta }^2\mu }{X^2}{y}_5^2-\frac{2E_p{\delta }^2}{X}{y}_5{y}_6+\frac{2E_p{\delta }^2}{X^3}{y}_4{y}_5-\frac{12E_p{\delta }^2}{X}{y}_4^3{y}_5\end{array}$

$\begin{array}{l}-\frac{3E_p{\delta }^2\mu }{X^4}{y}_4^2+\frac{12\mu }{X}{y}_4{y}_5+\frac{6E_p{\delta }^2\mu }{X^3}{y}_4{y}_5+\frac{12E_p{\delta }^2\mu }{X}{y}_4^2{\phi }_1-\frac{12E_p{\delta }^2\mu }{X^2}{y}_2{y}_4^2\\ -\frac{24E_p{\delta }^2}{X}{y}_2{y}_4{y}_5+\frac{24E_p{\delta }^2\mu }{X}{y}_2{\phi }_1-\frac{24E_p{\delta }^2\mu }{X^2}{y}_1{\phi }_1+\frac{72E_p{\delta }^2\mu }{X^3}{y}_1{y}_2+\frac{2\mu E_p{\delta }^2}{X}{y}_5{y}_6\\ +\frac{12E_p{\delta }^2\mu }{X}{y}_4^3{y}_5-\frac{2E_p{\delta }^2\mu }{X^2}{y}_4{y}_6+\frac{E_p{\delta }^2}{X^2}{y}_5^2-\frac{3E_p{\delta }^2}{X^4}{y}_4^2+\frac{3E_p{\delta }^2}{X^2}{y}_4^4-12E_p{\delta }^2{y}_4^3{y}_6\\ -36E_p{\delta }^2{y}_4^2{y}_5^2+\frac{12E_p{\delta }^2}{X^2}{y}_2^2-24E_p{\delta }^2{y}_2{y}_5^2-\frac{36E_p{\delta }^2}{X^4}{y}_1^2]\cdot 2E_p{\delta }^2{y}_5\end{array}$

$\begin{array}{l}{s}_4=[-4E_p{\delta }^2{y}_6{\kappa }_1-E_p{\delta }^2{y}_5^3-4E_p{\delta }^2{y}_4{y}_5{y}_6+\frac{1}{X^2}{y}_5-\frac{1}{X^3}{y}_4-\frac{2}{X}{y}_6-\frac{2E_p{\delta }^2\mu }{X^2}{y}_4{\phi }_1\\ -\frac{4E_p{\delta }^2\mu }{X^2}{y}_2{y}_5+\frac{2E_p{\delta }^2}{X^3}{y}_2{y}_4+\frac{4E_p{\delta }^2\mu }{X^3}{y}_2{y}_4-\frac{4E_p{\delta }^2}{X^4}{y}_1{y}_4-\frac{4E_p{\delta }^2\mu }{X^4}{y}_1{y}_4\\ +\frac{2E_p{\delta }^2}{X^3}{y}_1{y}_5+\frac{4E_p{\delta }^2\mu }{X^3}{y}_1{y}_5-\frac{2E_p{\delta }^2\mu }{X^2}{y}_1{y}_6+\frac{2E_p{\delta }^2\mu }{X}{y}_5{\phi }_1-\frac{4E_p{\delta }^2}{X}{y}_5{\phi }_1\\ +\frac{2E_p{\delta }^2\mu }{X}{y}_2{y}_6-\frac{4E_p{\delta }^2}{X}{y}_2{y}_6+\frac{2E_p{\delta }^2\mu }{X}{y}_4{y}_5^2-\frac{3E_p{\delta }^2}{X}{y}_4{y}_5^2+\frac{E_p{\delta }^2\mu }{X}{y}_4^2{y}_6\end{array}$

$\begin{array}{l}-\frac{2E_p{\delta }^2}{X}{y}_4^2{y}_6-\frac{E_p{\delta }^2\mu }{X^3}{y}_4^3+\frac{3E_p{\delta }^2\mu }{X^2}{y}_4^2{y}_5+\frac{12}{X}{y}_2{y}_4+\frac{12\mu }{X}{y}_2{y}_4+\frac{12E_p{\delta }^2}{X}{y}_2^2{y}_4\\ +\frac{12E_p{\delta }^2}{X}{y}_2{y}_4^3+\frac{6}{X}{y}_4^3+\frac{3E_p{\delta }^2}{X}{y}_4^5-\frac{12E_p{\delta }^2\mu }{X^3}{y}_1^2{y}_4+12y_2{y}_5+12E_p{\delta }^2{y}_2^2{y}_5\\ +36E_p{\delta }^2{y}_2{y}_4^2{y}_5+18y_4^2{y}_5+15E_p{\delta }^2{y}_4^4{y}_5+\frac{12\mu }{X}{y}_1{y}_5+\frac{12E_p{\delta }^2\mu }{X^2}{y}_1^2{y}_5+12y_4{\phi }_1\\ +24E_p{\delta }^2{y}_2{y}_4{\phi }_1+12E_p{\delta }^2{y}_4^3{\phi }_1+\frac{24E_p{\delta }^2\mu }{X^2}{y}_1{y}_2{y}_4+Q]\cdot (12+24E_p{\delta }^2{y}_2+12E_p{\delta }^2{y}_4^2)\end{array}$

$s_5=12+48E_P{\delta }^2{y}_2+24E_p{\delta }^2{y}_4^2+48E_p^2{\delta }^4{y}_2^2+48E_p^2{\delta }^4{y}_2{y}_4^2+12E_p^2{\delta }^4{y}_4^4-4E_p^2{\delta }^4{y}_5^4$

相应的边界条件：

周边简支：在 $x=1$ 处：

$y_1=0,y_3=0,-y_5-2E_P{\delta }^2{y}_2{y}_5-E_P{\delta }^2{y}_4^2{y}_5-\mu y_4=0$

周边固支：在 $x=1$ 处：

$y_1=0,y_3=0,y_4=0$

中心对称条件：在 $x=0$ 处

4. 数值求解及结果分析

在本文对周边固支的大挠度问题的分析中，考虑厚径比 $\delta =0.03$，本构非线性参数 $E_P$ 不同时，图1给出了圆板中心挠度 $W_0$ 随载荷Q变化的平衡路径，从图中不同 $E_P$ 值的路径曲线可以发现， $E_P=0$ 时挠度随载荷Q的增加而非线性增加， $E_P\ne 0$ 时挠度–载荷关系同样呈现出非线性，而且 $E_P$ 取正值和负值时，分居在 $E_P=0$ 曲线两侧。当 $E_P>0$ 时曲线在下方，挠度变小。这是由于大挠度弯曲变形时，面内挠度的增加使得弹性模量整体上增大而呈现出圆板的抗弯曲变形的能力增加了。当 $E_P<0$ 时，挠度曲线在上方，反映出圆板的抗弯曲变形能力的下降。

图2(a)给出了在 $Q=50$， $\delta =0.03$，本构非线性参数 $E_P$ 不同时，周边固支圆板的大挠度的弯曲构形图，从图中可以看出 $E_P$ 不同时板的弯曲构型基本相同但板的弯曲变形随着 $E_P$ 的增大而减小。由于几何中面变形的约束，使得横截面上的内力和弯矩的共同作用使圆板的中性层位置也发生了变化。图2(b)给出了在 $Q=50$， $\delta =0.03$ 时周边固支圆板的中性层位置。图中，纵坐标为弯曲变形时中性层的位置。由于几何非线性的影响，中性层位置不会因 $E_P$ 的正负变化而反向，这说明大挠度问题中，薄膜力对弯曲变形的影响更大一些。在 $x=0.6$ 附近反变点处弯矩很小，中性层远离几何中面，较大时中性层更靠近几何中面。

Figure 1. Equilibrium path for large deflection problem of circular plate with peripheral fixed support

图1. 周边固支圆板的大挠度问题的平衡路径

(a) 平衡构形 (b) 中性层位置

Figure 2. Equilibrium configuration and neutral layer position of peripheral fixed circular plate

图2. 周边固支圆板的平衡构形和中性层位置

Figure 3. Equilibrium path for large deflection problem of circumferential simply supported circular plate

图3. 周边简支圆板的大挠度问题的平衡路径

图3给出了周边简支圆板中心挠度 $W_0$ 随载荷Q变化的平衡路径，从图中可以发现， $E_P$ 不同时圆板的弯曲构形基本相同，圆板的弯曲变形同样随着 $E_P$ 的增加而减小，挠度随载荷Q的增加而非线性增加，而且 $E_P$ 取正值和负值时，分居在 $E_P=0$ 曲线两侧。当 $E_P>0$ 时曲线在下方，挠度变小，反应了圆板的抗弯曲变形能力的增加。当 $E_P<0$ 时，曲线在上方，则反应出圆板的抗弯曲变形能力的下降。

图4给出了在 $Q=50$， $\delta =0.03$，本构非线性参数 $E_P$ 不同时，周边简支圆板大挠度弯曲构形图和中性层位置。从图中可以看出，圆板周边简支时情形与周边固支时类似：随着 $E_P$ 的增大，板的抗弯曲变形能力的增加，中性层的位置也会相应的变化。

(a) 平衡构形 (b) 中性层位置

Figure 4. Equilibrium configuration and neutral layer position of simply supported circular plate

图4. 周边简支圆板的平衡构形和中性层位置

4. 结论

本文基于非线性本构关系，研究了周边固支和简支边界圆板在横向均布载荷的作用下的弯曲问题。假设材料弹性模量是应变的线性函数，运用数值法求解了该支撑条件下的大挠度问题的数值结果。结果表明：

1) 本构关系的非线性参数对板的弯曲变形有影响。弯曲变形与 $E_P$ 的符号有关并且随着的增大而减小。

2) 由于几何中面轴向变形的约束，使得 $E_P\ne 0$ 时截面上产生面内薄膜力，板的中性层位置发生变化，且截面中性层随着 $E_P$ 的变化而变化。

3) 非线性本构关系下的弯曲问题受本构非线性参数的影响，应采用大挠度问题的基本方程进行问题的求解更为妥当。

基金项目

本文得到国家自然科学基金11472123的资助。