拉伸载荷作用下复合材料干涉连接应力分析方法

摘要

本发明提出一种拉伸载荷作用下复合材料干涉连接应力分析方法，该方法基于Lekhnitskii复势函数理论，在建立接头孔边变形与应力函数关系的基础上，通过孔边变形的求解和应力函数运算，得到孔边各处的应力分量。该方法可以在已知材料属性、外载荷大小和干涉量的条件下求解接头孔边的应力状态，使用方便、计算准确，可以为结构强度的分析和工艺参数的分析提供有利的理论支持；弥补现有方法在计算时间和使用方便性上的缺陷，为结构的强度分析和损伤预测提供依据，为复合材料结构干涉配合连接的设计提供参考。

说明书

技术领域

本发明涉及复合材料连接应力分析技术领域，具体为一种拉伸载荷作用下复合材 料干涉连接应力分析方法。

背景技术

复合材料因其比强度高、比刚度高、可设计性强等优点，在航空、航天和汽车等工 业领域得到了越来越多的应用。复合材料在工程中的使用不可避免的导致连接结构的存 在。对于承载较大的结构而言，机械连接是最常用的连接形式。但由于机械连接需要预制连 接孔，造成材料不连续和应力集中，机械连接接头往往是结构中的最薄弱的区域，相关统计 结果表明60％～80％的结构失效发生在接头部位。为提高复合材料机械连接接头的强度和 疲劳寿命，美国麦道公司将已在金属结构中得到应用的干涉配合技术使用在复合材料结构 上，试验研究表明干涉配合对复合材料接头同样有效。但由于工艺要求较高且相关工艺参 数的选取较为复杂，目前该技术的应用尚处于发展状态。

拉伸载荷是连接结构承受的最主要外载形式之一。对于干涉配合的复合材料接头 而言，由于存在干涉残余应力并受到复合材料各向异性的影响，其连接孔周围的应力分布 沿角度变化剧烈，应力集中严重，极易成为损伤的源头，引起结构的破坏。应力水平的计算 是分析结构强度和寿命的基础，对工艺参数的选择具有重要作用。快速、准确地分析干涉连 接在外力作用下的应力状态，可以方便干涉量等工艺参数的选择，对于干涉配合技术的应 用具有十分重要的意义。

目前已有学者围绕复合材料连接接头的应力状态进行了研究，主要分为两类：间 隙配合接头在拉伸载荷作用下应力状态研究(AlukoOandWhitworthHA.“Analysisof stressdistributionaroundpinloadedholesinorthotropicplates”.Composite Structures.2008；86:308-13)；干涉连接接头周围干涉残余应力的研究(KimS,HeB,Shim C,andKimD,“Anexperimentalandnumericalstudyontheinterference-fitpin installationprocessforcross-plyglassfiberreinforcedplastics(GFRP)”, Composites:PartB.201354:153-162)。孔周应力状态的分析目前主要采用有限元数值模 拟的方式(张俊琪，刘龙权，陈昆昆，汪海.“干涉配合对复合材料机械连接结构承载能力的 影响”.上海交通大学学报,2013,47(11):1795-1806)。有限元数值模拟按照研究对象的几 何尺寸、材料属性和受力状态建立模型进行分析，尽管具有较高的精度但需要较长的计算 时间。而且由于模型与研究对象存在直接的对应关系，当研究对象的尺寸、材料、受力状态 发生变化时，需要重复建模和分析过程，因而对于大量不同工艺参数及材料属性的对象，其 所耗费的计算成本成倍增加。

由于承载能力强、便于拆卸等优点，飞机结构中大量使用机械连接接头。为了满足 寿命、强度、密封性等要求，干涉配合在复合材料机械连接接头中得到应用。对于飞机结构 件而言，外载作用下的应力状态是分析结构质量、安全性的关键因素，直接影响飞机产品的 安全可靠。目前，针对拉伸载荷作用下复合材料干涉配合接头的孔边应力状态主要采用有 限元数值模拟的方式，这种方法存在以下问题：1)计算时间代价很高。为保证精度，有限元 模型需要划分一定规模的网格单元，因而其计算过程较长。此外，提高计算精度需要成倍增 加计算规模，导致计算时间的进一步增加；2)模型具有一定的特异性。有限元模型针对特定 材料、特定尺寸的结构建立，当分析对象变化时，需要重复建模及计算过程，而该过程耗时 很长，使得其不适于大量对象的分析。

发明内容

为解决现有技术存在的问题，本发明提出了一种拉伸载荷作用下复合材料干涉连 接应力分析方法，该方法基于Lekhnitskii复势函数理论，在建立接头孔边变形与应力函数 关系的基础上，通过孔边变形的求解和应力函数运算，得到孔边各处的应力分量。该方法可 以在已知材料属性、外载荷大小和干涉量的条件下求解接头孔边的应力状态，使用方便、计 算准确，可以为结构强度的分析和工艺参数的分析提供有利的理论支持；弥补现有方法在 计算时间和使用方便性上的缺陷，为结构的强度分析和损伤预测提供依据，为复合材料结 构干涉配合连接的设计提供参考。

本发明的技术方案为：

所述一种拉伸载荷作用下复合材料干涉连接应力分析方法，其特征在于：复合材 料板上某一点在笛卡尔坐标系下的应力分量σ_x,σ_y,σ_xy通过下面公式表示，所述笛卡尔坐标 系为在复合材料板平面上，建立的以干涉连接孔中心为原点，以拉伸载荷方向为x轴正方向 的笛卡尔坐标系；复合材料板上该点在笛卡尔坐标系下的坐标为(x,y)；

$\left(\begin{array}{c}{\sigma }_x=2\mathrm{Re}\left[{\mu }_1^2{\Phi }_1^{\prime }\right(z_1)+{\mu }_2^2{\Phi }_2^{\prime }(z_2\left)\right]\\ {\sigma }_y=2\mathrm{Re}\left[{\Phi }_1^{\prime }\right(z_1)+{\Phi }_2^{\prime }(z_2\left)\right]\\ {\tau }_{xy}=-2\mathrm{Re}\left[{\mu }_1{\Phi }_1^{\prime }\right(z_1)+{\mu }_2{\Phi }_2^{\prime }(z_2\left)\right]\end{array}\right)$

其中Re表示取实部，μ₁、μ₂为与复合材料属性相关的复参数，Φ′₁(z₁)为应力函数 Φ₁(z₁)的导数，Φ′₂(z₂)为应力函数Φ₂(z₂)的导数；

μ₁、μ₂通过公式

$\left(\begin{array}{c}{\mu }_1{\mu }_2+k=0\\ i({\mu }_1+{\mu }_2)+n=0\end{array}\right)$

得到，其中i为虚数单位，k和n为与复合材料参数相关的系数，k＝(E₁/E₂)^1/2，n＝[2 (k-v₁₂)+E₁/G₁₂]^1/2，E₁、E₂为复合材料弹性模量，G₁₂为复合材料剪切模量；

应力函数Φ₁(z₁)和Φ₂(z₂)为：

应力函数Φ₁(z₁)和Φ₂(z₂)的导数为：

$\left(\begin{array}{c}{\Phi }_1^{\prime }\left(z_1\right)=[A-\frac{\Delta }{2D}\frac{q_2-{\mathrm{ip}}_2}{{\zeta }_1}-\frac{1}{D}(U_1{q}_2-{\mathrm{iV}}_1{p}_2)\frac{1}{{\zeta }_1^2}-\frac{2}{D}(U_2{q}_2-{\mathrm{iV}}_2{p}_2)\frac{1}{{\zeta }_1^4}]\frac{1}{\sqrt{z_1^2-{\mu }_1{r}^2-r^2}}\\ {\Phi }_2^{\prime }\left(z_2\right)=[B+\frac{\Delta }{2D}\frac{q_1-{\mathrm{ip}}_1}{{\zeta }_2}+\frac{1}{D}(U_1{q}_1-{\mathrm{iV}}_1{p}_1)\frac{1}{{\zeta }_2^2}+\frac{2}{D}(U_2{q}_2-{\mathrm{iV}}_2{p}_1)\frac{1}{{\zeta }_2^4}]\frac{1}{\sqrt{z_2^2-{\mu }_2{r}^2-r^2}}\end{array}\right)$

其中A和B是与拉伸载荷以及复合材料参数相关的系数，

$\left(\begin{array}{c}A=\frac{P}{\pi i}\frac{{\mu }_1{\bar{\mu }}_1+{\mu }_1{\mu }_2+{\mu }_1{\bar{\mu }}_2-(a_{12}/a_{22}){\mu }_1{\mu }_2{\bar{\mu }}_1{\bar{\mu }}_2}{({\mu }_1-{\bar{\mu }}_1)({\mu }_1-{\mu }_2)({\mu }_1-{\bar{\mu }}_2)}\\ B=\frac{P}{\pi i}\frac{{\mu }_2{\bar{\mu }}_2+{\mu }_2{\mu }_1+{\mu }_2{\bar{\mu }}_1-(a_{12}/a_{22}){\mu }_1{\mu }_2{\bar{\mu }}_1{\bar{\mu }}_2}{({\mu }_2-{\bar{\mu }}_2)({\mu }_2-{\mu }_1)({\mu }_2-{\bar{\mu }}_1)}\end{array}\right)$

P为拉伸载荷，和对应为μ₁和μ₂的共轭，Δ＝R-r表示干涉量，R为连接螺栓半 径，r为复合材料板干涉连接孔初始直径，D＝(μ₁-μ₂)g/E₁，g为与复合材料参数相关的系数， g＝(1-v₁₂v₂₁)/E₂+k/G₁₂，v₁₂、v₂₁为复合材料的泊松系数；

系数U₁，U₂，V₁，V₂通过以下公式得到：

$\left(\begin{array}{c}{U}_1=\frac{c-1}{2c}{u}_0\\ U_2=\frac{c+1}{2c}{u}_0=V_2\\ V_1=(\frac{c-1}{2c}+\frac{c+1}{c})u_0\end{array}\right)$

其中u₀为与拉伸载荷大小相关的位移量，c为与复合材料参数相关的系数，

$u_0=\frac{Pg}{4\pi }\frac{10(k-v_{12}+nk+2n)-11n(1-k)}{5n^{2}k(1-k)+[5(k-v_{12}+nk+2n)-3n(1-k)](k-v_{12}+nk)}$

c＝-(10(k-v₁₂+nk+2n)-11n(1-k))/(10(k-v₁₂+nk+2n)-n(1-k))

ζ_j为与坐标有关的参数，j＝1,2：

${\zeta }_j=\frac{z_j+\sqrt{z_j^2-{\mu }_j{r}^2-r^2}}{r-{\mathrm{i\mu }}_{j}r},z_j=x+{\mu }_{j}y$

p₁,p₂,q₁,q₂为与材料参数有关的系数，

$\left(\begin{array}{c}{p}_1=a_{11}{\mu }_1^2+a_{12}\\ p_2=a_{12}{\mu }_2^2+a_{12}\\ q_1=a_{12}{\mu }_1+\frac{a_{22}}{{\mu }_1}\\ q_2=a_{12}{\mu }_2+\frac{a_{22}}{{\mu }_2}\end{array}\right)$

其中a₁₁、a₁₂、a₂₂为复合材料板的柔度系数。

进一步的优选方案，所述一种拉伸载荷作用下复合材料干涉连接应力分析方法， 其特征在于：复合材料板上某一点在极坐标系下的应力分量σ_r,σ_θ,τ_rθ通过下面公式表示，

$\left(\begin{array}{c}{\sigma }_r=2\mathrm{Re}\left[{\left(\sin \theta -{\mu }_1\cos \theta \right)}^2{\Phi }_1^{\prime }\left(z_1\right)+{\left(\sin \theta -{\mu }_2\cos \theta \right)}^2{\Phi }_2^{\prime }\left(z_2\right)\right]\\ {\sigma }_{\theta }=2\mathrm{Re}\left[{\left({\mu }_1\sin \theta +\cos \theta \right)}^2{\Phi }_1^{\prime }\left(z_1\right)+{\left({\mu }_2\sin \theta +\cos \theta \right)}^2{\Phi }_2^{\prime }\left(z_2\right)\right]\\ {\tau }_{r\theta }=2\mathrm{Re}\left[\left(\sin \theta -{\mu }_1\cos \theta \right)\left(\cos \theta +{\mu }_1\sin \theta \right){\Phi }_1^{\prime }\left(z_1\right)+\left(\sin \theta -{\mu }_2\cos \theta \right)\left(\cos \theta +{\mu }_2\sin \theta \right){\Phi }_2^{\prime }\left(z_2\right)\right]\end{array}\right)$

所述极坐标系为在复合材料板平面上，建立的以干涉连接孔中心为原点，以拉伸 载荷方向为0°方向建立的极坐标系，θ表示点在极坐标系中的角度。

有益效果

本发明具有的优点如下：

(1)本发明建立了拉伸载荷作用下的复合材料干涉配合接头应力分布的计算模 型，可以在Matlab软件中编程求解，具有精度高、速度快、计算代价低、使用方便等优点。

(2)该方法具有一定的通用性，对于不同材料属性、不同干涉量的接头只需修改材 料及结构参数，不需重复建模及编程，可以降低分析工作的耗时。

(3)该方法可以实现整个复合材料板平面内任何点位的应力计算，可以为结构设 计及强度预测提供支持。

附图说明

本发明的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得 明显和容易理解，其中：

图1：结构受力状态及主要几何尺寸；

1表示金属螺栓，2表示复合材料板；

图2：螺栓孔边区域变形位移分析；

图3：拉伸载荷作用下螺栓孔边应力分布结果。

具体实施方式

下面详细描述本发明的实施例，所述实施例的示例在附图中示出，其中自始至终 相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附 图描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。

本发明的目的在于提出一种拉伸载荷作用下复合材料干涉连接应力分析方法，弥 补现有方法在计算时间和使用方便性上的缺陷，为结构的强度分析和损伤预测提供依据， 为复合材料结构干涉配合连接的设计提供参考。该方法计算成本低、计算量小并具有较高 的精度，可以满足工程应用中的需求。

本发明方法涉及到的主要部件为螺栓(1)、复合材料板(2)，其中螺栓半径为R、复 合材料板螺栓孔初始直径为r(r<R)、复合材料板厚为h、复合材料宽度为W、螺栓孔边距为E。 螺栓通过冷挤压方式安装到复合材料螺栓孔中。

申请人基于上述研究目的，提出本发明的研究思路如下：

(1)根据拉伸外载荷与孔周区域变形的关系，分析拉伸载荷引起的孔周变形位移。 由于拉伸载荷受力沿板的长度方向对称，因而仅需分析模型的一半。根据变形特点，可以将 孔周按角度分为5个区域：区域Ⅰ为连接件与复合材料板的纯挤压区域，区域Ⅱ为连接件与 复合材料板相互挤压并存在滑动的区域，区域Ⅲ为连接件与复合材料板间不存在法向挤压 的区域，区域Ⅳ为连接件与复合材料板相互分离并存在滑动的区域，区域Ⅴ为连接件与复 合材料板相互分离的区域。忽略螺栓的变形，以孔中心为原点，拉伸载荷方向为0°方向建立 极坐标系，用u、v表示沿0°和90°的位移分量，各区域的范围及位移关系如下：

区域Ⅰθ＝0u＝u₀andv＝0

区域Ⅱ0<θ<π/2(u₀-u)cosθ＝vsinθ

区域Ⅲθ＝π/2u＝u₀/candv＝0

区域Ⅳπ/2<θ<π(u₀-u)cosθ＝-vsinθandσ_r≤0

区域Ⅴθ＝πu＝u₀andσ_r≤0

其中σ_r为某一点的径向应力，θ表示孔上一点在坐标系中的角度(逆时针旋转为 正)，u₀为与外载荷大小相关的位移量，其数值可由公式(1)给出

$u_0=\frac{Pg}{4\pi }\frac{10(k-v_{12}+nk+2n)-11n(1-k)}{5n^{2}k(1-k)+[5(k-v_{12}+nk+2n)-3n(1-k)](k-v_{12}+nk)}---\left(1\right)$

P为外载荷，k,n,g,c均为与材料参数相关的系数，其计算公式为

k＝(E₁/E₂)^1/2,

n＝[2(k-v₁₂)+E₁/G₁₂]^1/2,

g＝(1-v₁₂v₂₁)/E₂+k/G₁₂,

c＝-(10(k-v₁₂+nk+2n)-11n(1-k))/(10(k-v₁₂+nk+2n)-n(1-k))。

计算公式中v₁₂，v₂₁为复合材料的泊松系数，E₁，E₂为弹性模量，G₁₂为剪切模量。最终 可以将位移分量u和v表示为公式(2)的形式

$\left(\begin{array}{c}u=U_{1}cos2\theta +U_{2}cos4\theta \\ v=V_1\sin 2\theta +V_2\sin 4\theta \end{array}\right)---\left(2\right)$

系数U₁，U₂，V₁，V₂用于确定应力函数。在(2)式带入区域I、Ⅱ、Ⅲ的位移关系得如下 关系式

$\left(\begin{array}{c}{U}_1=\frac{c-1}{2c}{u}_0\\ U_2=\frac{c+1}{2c}{u}_0=V_2\\ V_1=(\frac{c-1}{2c}+\frac{c+1}{c})u_0\end{array}\right)---\left(3\right)$

即U₁，U₂，V₁，V₂均可由u₀和c的值得到。

(2)对于受到拉伸载荷作用的干涉配合接头而言，其应力状态可以视为干涉应力 与拉伸应力的叠加。在板平面上，以孔中心为原点，拉伸载荷方向为x轴正方向建立笛卡尔 坐标系。依据Lekhnitskii复势函数理论，根据变形位移与应力状态的关系，其应力函数Φ 表达式为

$\left(\begin{array}{c}{\Phi }_1\left(z_1\right)={\mathrm{Aln\zeta }}_1+\frac{\Delta }{2D}\frac{q_2-{\mathrm{ip}}_2}{{\zeta }_1}+\frac{1}{2D}(U_1{q}_2-{\mathrm{iV}}_1{p}_2)\frac{1}{{\zeta }_1^2}+\frac{1}{2D}(U_2{q}_2-{\mathrm{iV}}_2{p}_2)\frac{1}{{\zeta }_1^4}\\ {\Phi }_2\left(z_2\right)={\mathrm{Bln\zeta }}_2-\frac{\Delta }{2D}\frac{q_1-{\mathrm{ip}}_1}{{\zeta }_2}-\frac{1}{2D}(U_1{q}_1+{\mathrm{iV}}_1{p}_1)\frac{1}{{\zeta }_2^2}-\frac{1}{2D}(U_2{q}_1+{\mathrm{iV}}_2{p}_1)\frac{1}{{\zeta }_2^4}\end{array}\right)---\left(4\right)$

式中Δ＝R-r表示干涉量，D＝(μ₁-μ₂)g/E₁。μ₁，μ₂为与材料属性相关的复参数，其值 由式(5)求得。

$\left(\begin{array}{c}{\mu }_1{\mu }_2+k=0\\ i({\mu }_1+{\mu }_2)+n=0\end{array}\right)---\left(5\right)$

ζ_j(j＝1,2)为与坐标有关的参数。对复合材料板上任意一点(x,y)而言，其坐标的 复参数表达形式为z_j＝x+μ_jy，对该表达式进行保角变化得ζ_j，如公式(6)

${\zeta }_j=\frac{z_j+\sqrt{z_j^2-{\mu }_j{r}^2-r^2}}{r-{\mathrm{i\mu }}_{j}r}---\left(6\right)$

p₁,p₂,q₁,q₂为与材料参数有关的系数，其值可通过方程组(7)求解得到。a_mn(m,n＝ 1,2,6)为复合材料板的柔度系数。

$\left(\begin{array}{c}{p}_1=a_{11}{\mu }_1^2+a_{12}\\ p_2=a_{12}{\mu }_2^2+a_{12}\\ q_1=a_{12}{\mu }_1+\frac{a_{22}}{{\mu }_1}\\ q_2=a_{12}{\mu }_2+\frac{a_{22}}{{\mu }_2}\end{array}\right)---\left(7\right)$

A和B是与外载荷和材料参数相关的系数，其表达式如公式(8),i为虚数单位。

$\left(\begin{array}{c}A=\frac{P}{\pi i}\frac{{\mu }_1{\bar{\mu }}_1+{\mu }_1{\mu }_2+{\mu }_1{\bar{\mu }}_2-(a_{12}/a_{22}){\mu }_1{\mu }_2{\bar{\mu }}_1{\bar{\mu }}_2}{({\mu }_1-{\bar{\mu }}_1)({\mu }_1-{\mu }_2)({\mu }_1-{\bar{\mu }}_2)}\\ B=\frac{P}{\pi i}\frac{{\mu }_2{\bar{\mu }}_2+{\mu }_2{\mu }_1+{\mu }_2{\bar{\mu }}_1-(a_{12}/a_{22}){\mu }_1{\mu }_2{\bar{\mu }}_1{\bar{\mu }}_2}{({\mu }_2-{\bar{\mu }}_2)({\mu }_2-{\mu }_1)({\mu }_2-{\bar{\mu }}_1)}\end{array}\right)---\left(8\right)$

(3)由于应力分量是由应力函数的导数表达的，对(3)式求导并进一步化简得到

$\left(\begin{array}{c}{\Phi }_1^{\prime }\left(z_1\right)=[A-\frac{\Delta }{2D}\frac{q_2-{\mathrm{ip}}_2}{{\zeta }_1}-\frac{1}{D}(U_1{q}_2-{\mathrm{iV}}_1{p}_2)\frac{1}{{\zeta }_1^2}-\frac{2}{D}(U_2{q}_2-{\mathrm{iV}}_2{p}_2)\frac{1}{{\zeta }_1^4}]\frac{1}{\sqrt{z_1^2-{\mu }_1{r}^2-r^2}}\\ {\Phi }_2^{\prime }\left(z_2\right)=[B+\frac{\Delta }{2D}\frac{q_1-{\mathrm{ip}}_1}{{\zeta }_2}+\frac{1}{D}(U_1{q}_1-{\mathrm{iV}}_1{p}_1)\frac{1}{{\zeta }_2^2}+\frac{2}{D}(U_2{q}_2-{\mathrm{iV}}_2{p}_1)\frac{1}{{\zeta }_2^4}]\frac{1}{\sqrt{z_2^2-{\mu }_2{r}^2-r^2}}\end{array}\right)---\left(9\right)$

(4)在笛卡尔坐标系下的应力分量σ_x,σ_y,σ_xy可由应力函数的导数公式(9)求得，表 达式如(10)所示。考虑到接头孔的形状特点，应力分量采用极坐标的表达意义更为清晰，便 于理解。借助转轴公式可以得到极坐标系下的应力分量σ_r,σ_θ,τ_rθ表达式(11)

$\left(\begin{array}{c}{\sigma }_x=2\mathrm{Re}\left[{\mu }_1^2{\Phi }_1^{\prime }\right(z_1)+{\mu }_2^2{\Phi }_2^{\prime }(z_2\left)\right]\\ {\sigma }_y=2\mathrm{Re}\left[{\Phi }_1^{\prime }\right(z_1)+{\Phi }_2^{\prime }(z_2\left)\right]\\ {\tau }_{xy}=-2\mathrm{Re}\left[{\mu }_1{\Phi }_1^{\prime }\right(z_1)+{\mu }_2{\Phi }_2^{\prime }(z_2\left)\right]\end{array}\right)---\left(10\right)$

$\left(\begin{array}{c}{\sigma }_r=2\mathrm{Re}\left[{\left(\sin \theta -{\mu }_1\cos \theta \right)}^2{\Phi }_1^{\prime }\left(z_1\right)+{\left(\sin \theta -{\mu }_2\cos \theta \right)}^2{\Phi }_2^{\prime }\left(z_2\right)\right]\\ {\sigma }_{\theta }=2\mathrm{Re}\left[{\left({\mu }_1\sin \theta +\cos \theta \right)}^2{\Phi }_1^{\prime }\left(z_1\right)+{\left({\mu }_2\sin \theta +\cos \theta \right)}^2{\Phi }_2^{\prime }\left(z_2\right)\right]\\ {\tau }_{r\theta }=2\mathrm{Re}\left[\left(\sin \theta -{\mu }_1\cos \theta \right)\left(\cos \theta +{\mu }_1\sin \theta \right){\Phi }_1^{\prime }\left(z_1\right)+\left(\sin \theta -{\mu }_2\cos \theta \right)\left(\cos \theta +{\mu }_2\sin \theta \right){\Phi }_2^{\prime }\left(z_2\right)\right]\end{array}\right)---\left(11\right)$

。

基于上述研究过程，本发明的方法为复合材料板上某一点在笛卡尔坐标系下的应 力分量σ_x,σ_y,σ_xy通过下面公式表示，所述笛卡尔坐标系为在复合材料板平面上，建立的以 干涉连接孔中心为原点，以拉伸载荷方向为x轴正方向的笛卡尔坐标系；复合材料板上该点 在笛卡尔坐标系下的坐标为(x,y)；

$\left(\begin{array}{c}{\sigma }_x=2\mathrm{Re}\left[{\mu }_1^2{\Phi }_1^{\prime }\right(z_1)+{\mu }_2^2{\Phi }_2^{\prime }(z_2\left)\right]\\ {\sigma }_y=2\mathrm{Re}\left[{\Phi }_1^{\prime }\right(z_1)+{\Phi }_2^{\prime }(z_2\left)\right]\\ {\tau }_{xy}=-2\mathrm{Re}\left[{\mu }_1{\Phi }_1^{\prime }\right(z_1)+{\mu }_2{\Phi }_2^{\prime }(z_2\left)\right]\end{array}\right)$

其中Re表示取实部，μ₁、μ₂为与复合材料属性相关的复参数，Φ′₁(z₁)为应力函数 Φ₁(z₁)的导数，Φ′₂(z₂)为应力函数Φ₂(z₂)的导数；

μ₁、μ₂通过公式

$\left(\begin{array}{c}{\mu }_1{\mu }_2+k=0\\ i({\mu }_1+{\mu }_2)+n=0\end{array}\right)$

得到，其中i为虚数单位，k和n为与复合材料参数相关的系数，k＝(E₁/E₂)^1/2，n＝[2 (k-v₁₂)+E₁/G₁₂]^1/2，E₁、E₂为复合材料弹性模量，G₁₂为复合材料剪切模量；

应力函数Φ₁(z₁)和Φ₂(z₂)为：

应力函数Φ₁(z₁)和Φ₂(z₂)的导数为：

$\left(\begin{array}{c}{\Phi }_1^{\prime }\left(z_1\right)=[A-\frac{\Delta }{2D}\frac{q_2-{\mathrm{ip}}_2}{{\zeta }_1}-\frac{1}{D}(U_1{q}_2-{\mathrm{iV}}_1{p}_2)\frac{1}{{\zeta }_1^2}-\frac{2}{D}(U_2{q}_2-{\mathrm{iV}}_2{p}_2)\frac{1}{{\zeta }_1^4}]\frac{1}{\sqrt{z_1^2-{\mu }_1{r}^2-r^2}}\\ {\Phi }_2^{\prime }\left(z_2\right)=[B+\frac{\Delta }{2D}\frac{q_1-{\mathrm{ip}}_1}{{\zeta }_2}+\frac{1}{D}(U_1{q}_1-{\mathrm{iV}}_1{p}_1)\frac{1}{{\zeta }_2^2}+\frac{2}{D}(U_2{q}_2-{\mathrm{iV}}_2{p}_1)\frac{1}{{\zeta }_2^4}]\frac{1}{\sqrt{z_2^2-{\mu }_2{r}^2-r^2}}\end{array}\right)$

其中A和B是与拉伸载荷以及复合材料参数相关的系数，

$\left(\begin{array}{c}A=\frac{P}{\pi i}\frac{{\mu }_1{\bar{\mu }}_1+{\mu }_1{\mu }_2+{\mu }_1{\bar{\mu }}_2-(a_{12}/a_{22}){\mu }_1{\mu }_2{\bar{\mu }}_1{\bar{\mu }}_2}{({\mu }_1-{\bar{\mu }}_1)({\mu }_1-{\mu }_2)({\mu }_1-{\bar{\mu }}_2)}\\ B=\frac{P}{\pi i}\frac{{\mu }_2{\bar{\mu }}_2+{\mu }_2{\mu }_1+{\mu }_2{\bar{\mu }}_1-(a_{12}/a_{22}){\mu }_1{\mu }_2{\bar{\mu }}_1{\bar{\mu }}_2}{({\mu }_2-{\bar{\mu }}_2)({\mu }_2-{\mu }_1)({\mu }_2-{\bar{\mu }}_1)}\end{array}\right)$

P为拉伸载荷，和对应为μ₁和μ₂的共轭，Δ＝R-r表示干涉量，R为连接螺栓半 径，r为复合材料板干涉连接孔初始直径，D＝(μ₁-μ₂)g/E₁，g为与复合材料参数相关的系数， g＝(1-v₁₂v₂₁)/E₂+k/G₁₂，v₁₂、v₂₁为复合材料的泊松系数；

系数U₁，U₂，V₁，V₂通过以下公式得到：

$\left(\begin{array}{c}{U}_1=\frac{c-1}{2c}{u}_0\\ U_2=\frac{c+1}{2c}{u}_0=V_2\\ V_1=(\frac{c-1}{2c}+\frac{c+1}{c})u_0\end{array}\right)$

其中u₀为与拉伸载荷大小相关的位移量，c为与复合材料参数相关的系数，

$u_0=\frac{Pg}{4\pi }\frac{10(k-v_{12}+nk+2n)-11n(1-k)}{5n^{2}k(1-k)+[5(k-v_{12}+nk+2n)-3n(1-k)](k-v_{12}+nk)}$

c＝-(10(k-v₁₂+nk+2n)-11n(1-k))/(10(k-v₁₂+nk+2n)-n(1-k))

ζ_j为与坐标有关的参数，j＝1,2：

${\zeta }_j=\frac{z_j+\sqrt{z_j^2-{\mu }_j{r}^2-r^2}}{r-{\mathrm{i\mu }}_{j}r},z_j=x+{\mu }_{j}y$

p₁,p₂,q₁,q₂为与材料参数有关的系数，

$\left(\begin{array}{c}{p}_1=a_{11}{\mu }_1^2+a_{12}\\ p_2=a_{12}{\mu }_2^2+a_{12}\\ q_1=a_{12}{\mu }_1+\frac{a_{22}}{{\mu }_1}\\ q_2=a_{12}{\mu }_2+\frac{a_{22}}{{\mu }_2}\end{array}\right)$

其中a₁₁、a₁₂、a₂₂为复合材料板的柔度系数。

进一步的，复合材料板上某一点在极坐标系下的应力分量σ_r,σ_θ,τ_rθ通过下面公式 表示，

$\left(\begin{array}{c}{\sigma }_r=2\mathrm{Re}\left[{\left(\sin \theta -{\mu }_1\cos \theta \right)}^2{\Phi }_1^{\prime }\left(z_1\right)+{\left(\sin \theta -{\mu }_2\cos \theta \right)}^2{\Phi }_2^{\prime }\left(z_2\right)\right]\\ {\sigma }_{\theta }=2\mathrm{Re}\left[{\left({\mu }_1\sin \theta +\cos \theta \right)}^2{\Phi }_1^{\prime }\left(z_1\right)+{\left({\mu }_2\sin \theta +\cos \theta \right)}^2{\Phi }_2^{\prime }\left(z_2\right)\right]\\ {\tau }_{r\theta }=2\mathrm{Re}\left[\left(\sin \theta -{\mu }_1\cos \theta \right)\left(\cos \theta +{\mu }_1\sin \theta \right){\Phi }_1^{\prime }\left(z_1\right)+\left(\sin \theta -{\mu }_2\cos \theta \right)\left(\cos \theta +{\mu }_2\sin \theta \right){\Phi }_2^{\prime }\left(z_2\right)\right]\end{array}\right)$

所述极坐标系为在复合材料板平面上，建立的以干涉连接孔中心为原点，以拉伸 载荷方向为0°方向建立的极坐标系，θ表示点在极坐标系中的角度。

本实施例中，复合材料平板通过钛合金高锁螺栓连接，连接采用干涉配合，结构的 材料属性及几何参数如表1、表2、表3所示，所受外载为匀布拉伸载荷P＝500Mpa。

表1干涉配合接头几何尺寸(mm)

注：E为螺栓孔中心到复合材料构件一端的最短距离，W为复合材料板的宽度，h为 复合材料板的厚度，r为复合材料板预制连接孔的半径，R为螺栓半径，如图1所示。

表2螺栓材料及其性能参数

注：剪切模量G可以通过弹性模量E和泊松比v求得。

表3复合材料板材料及其性能参数

注：E₁,E₂,E₃为复合材料的弹性模量，ν₁₂,ν₁₃,ν₂₃为复合材料的泊松比，G₁₂,G₁₃,G₂₃为复合材料的剪切模量。

拉伸载荷作用下的接头孔边应力计算具体实施步骤如下：

(1)将材料弹性模量、剪切模量、泊松比带入系数计算公式得

k＝3.8357,g＝0.0011,n＝6.2475

将系数值、材料属性和外载荷数值带入位移计算公式得到

u₀＝0.00525618590918473

c＝-1.424378995998533

U₁＝0.0044731727，U₂＝0.0007830131

V₁＝0.0060391990，V₂＝0.0007830131

孔周区域位移分量u，v为

$\left(\begin{array}{c}u=0.0044731727\cos 2\theta +0.0007830131\cos 4\theta \\ v=0.0060391990\sin 2\theta +0.0007830131\sin 4\theta \end{array}\right).$

(2)将材料参数带入复参数μ₁，μ₂的计算公式得

μ₁＝5.55730240104i

μ₂＝0.69021012986i

将材料参数和外载荷带入系数A、B的计算公式得

A＝14.6944803221

B＝-3.0196589997

将材料常数带入系数D，p₁,p₂,q₁,q₂的计算公式得

D＝0.0000000408025506045i

p₁＝-0.0002437782

p₂＝-0.0000062217

q₁＝(-0.0000345764)i

q₂＝(-0.0001682581)i

(3)将系数带入应力函数的导数表达式得

$\left(\begin{array}{c}{\Phi }_1^{\prime }\left(z_1\right)=[14.6944803221-\frac{0.0006\times \left(\right(-0.0001682581)\times i+0.0000062217\times i)}{2\times 0.0000000408025506045\times i}(\cos \theta -i\sin \theta )-\\ \frac{-0.0044731727\times 0.0001682581\times i+0.0060391990\times 0.0000062217\times i}{0.0000000408025506045\times i}(\cos 2\theta -i\sin 2\theta )-\\ \frac{2(-0.0044731727\times 0.0001682581\times i+0.0060391990\times 0.0000062217\times i)}{0.0000000408025506045\times i}(\cos 4\theta -i\sin 4\theta )]\\ \times \frac{1}{i\times 3\times (\sin \theta -5.55730240104\times i\times \cos \theta )}\\ {\Phi }_2^{\prime }\left(z_2\right)=[-3.0196589997+\frac{0.0006\times \left(\right(-0.0000345764)\times i+0.0002437782\times i)}{2\times 0.0000000408025506045\times i}(\cos \theta -i\sin \theta )-\\ \frac{-0.0044731727\times 0.0000345764\times i+0.0060391990\times 0.0002437782\times i}{0.0000000408025506045\times i}(\cos 2\theta -i\sin 2\theta )-\\ \frac{2(-0.0044731727\times 0.0000345764\times i+0.0060391990\times 0.0002437782\times i)}{0.0000000408025506045\times i}(\cos 4\theta -i\sin 4\theta )]\\ \times \frac{1}{i\times 3\times (\sin \theta -0.69021012986\times i\times \cos \theta )}\end{array}\right)$

(4)笛卡尔坐标系下的应力分量表达式如下：

$\left(\begin{array}{c}{\sigma }_x=2\mathrm{Re}\left[{\left(5.55730240104\times i\right)}^2{\Phi }_1^{\prime }\right(z_1)+{\left(0.69021012986\times i\right)}^2{\Phi }_2^{\prime }(z_2\left)\right]\\ {\sigma }_y=2\mathrm{Re}\left[{\Phi }_1^{\prime }\right(z_1)+{\Phi }_2^{\prime }(z_2\left)\right]\\ {\tau }_{xy}=-2\mathrm{Re}[5.55730240104\times i\times {\Phi }_1^{\prime }(z_1)+0.69021012986\times i\times {\Phi }_2^{\prime }(z_2\left)\right]\end{array}\right)$

极坐标系下的应力分量表达式如下

$\left(\begin{array}{c}{\sigma }_r=2\mathrm{Re}\left[{\left(\sin \theta -5.55730240104\times i\times \cos \theta \right)}^2{\Phi }_1^{\prime }\left(z_1\right)+{\left(\sin \theta -0.69021012986\times i\times \cos \theta \right)}^2{\Phi }_2^{\prime }\left(z_2\right)\right]\\ {\sigma }_{\theta }=2\mathrm{Re}\left[{\left(5.55730240104\times i\times \sin \theta +\cos \theta \right)}^2{\Phi }_1^{\prime }\left(z_1\right)+{\left(0.69021012986\times i\times \sin \theta +\cos \theta \right)}^2{\Phi }_2^{\prime }\left(z_2\right)\right]\\ {\tau }_{r\theta }=2\mathrm{Re}[\left(\sin \theta -5.55730240104\times i\times \cos \theta \right)\left(\cos \theta +5.55730240104\times i\times \sin \theta \right){\Phi }_1^{\prime }\left(z_1\right)+(\sin \theta \\ -0.69021012986\times i\times \cos \theta )\left(\cos \theta +0.69021012986\times i\times \sin \theta \right){\Phi }_2^{\prime }\left(z_2\right)]\end{array}\right)$

将应力函数导数结果带入应力分量表达式即可获得复合材料板孔周区域应力分 量的计算公式，带入所求位置的角度即可得到相应点位的应力状态。有限元模拟结果对比 如图3所示，由图可知计算得到的应力分布与有限元数值仿真结果一致，且得到的应力数值 较为精确，可以满足工程需要。

尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例 性的，不能理解为对本发明的限制，本领域的普通技术人员在不脱离本发明的原理和宗旨 的情况下在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。