三参数 S-N 曲线拟合及MATLAB程序

更多阅读：sppy.site

S-N 曲线：应力 S 与疲劳寿命 N 之间关系的曲线。

三参数 S-N 曲线方程
(S−Sf)m×N=C(1)\tag{1} (S-S_\mathrm{f})^m\times N=C (S−Sf)m×N=C(1)

线性拟合方法

将式(1)取对数，可得

mlg⁡(S−Sf)+lg⁡N=lg⁡C(2)\tag{2} m\lg(S-S_\mathrm{f})+\lg N=\lg C mlg(S−Sf)+lgN=lgC(2)
令 x=lg⁡(S−Sf)x=\lg(S-S_\mathrm{f})x=lg(S−Sf) 、y=lg⁡Ny=\lg Ny=lgN 、a=lg⁡Ca=\lg Ca=lgC 、b=−mb=-mb=−m ，可得
y=a+bx(3)\tag{3} y=a+bx y=a+bx(3)
对于一系列的应力值 SiS_iSi 、疲劳寿命值 NiN_iNi (i=1,2,3,⋯,ni=1,2,3,\cdots,ni=1,2,3,⋯,n)，由最小二乘法可得
{b=Lxy/Lxxa=yˉ−xˉb(4a)\tag{4a} \begin{cases} b=L_{xy}/L_{xx}\\[5pt] a=\bar{y}-\bar{x}b \end{cases} ⎩⎨⎧b=Lxy/Lxxa=yˉ−xˉb(4a)
线性相关系数 RRR 的平方为
R2=Lxy2LxxLyy(4b)\tag{4b} R^2=\frac{L_{xy}^2}{L_{xx}L_{yy}} R2=LxxLyyLxy2(4b)
其中
xi=lg⁡(Si−Sf)yi=lg⁡Nixˉ=1n∑i=1nxiyˉ=1n∑i=1nyiLxx=∑i=1nxi2−1n(∑i=1nxi)2Lyy=∑i=1nyi2−1n(∑i=1nyi)2Lxy=∑i=1nxiyi−1n(∑i=1nxi)(∑i=1nyi)x_i=\lg(S_i-S_\mathrm{f}) \qquad y_i=\lg N_i\\[8pt] \bar{x}=\frac{1}{n}\sum_{i=1}^nx_i \qquad \bar{y}=\frac{1}{n}\sum_{i=1}^ny_i\\[8pt] L_{xx}=\sum_{i=1}^nx_i^2-\frac{1}{n}\Big(\sum_{i=1}^nx_i\Big)^2\\[8pt] L_{yy}=\sum_{i=1}^ny_i^2-\frac{1}{n}\Big(\sum_{i=1}^ny_i\Big)^2\\[8pt] L_{xy}=\sum_{i=1}^nx_iy_i-\frac{1}{n}\Big(\sum_{i=1}^nx_i\Big)\Big(\sum_{i=1}^ny_i\Big) xi=lg(Si−Sf)yi=lgNixˉ=n1i=1∑nxiyˉ=n1i=1∑nyiLxx=i=1∑nxi2−n1(i=1∑nxi)2Lyy=i=1∑nyi2−n1(i=1∑nyi)2Lxy=i=1∑nxiyi−n1(i=1∑nxi)(i=1∑nyi)
由于 SfS_\mathrm{f}Sf 未知，故式(4)中的表达式均可看作关于 SfS_\mathrm{f}Sf 的函数。为了得到最佳的线性相关性，即线性相关系数 RRR 的绝对值最大，于是有
d(R2)dSf=0(5a)\tag{5a} \frac{\mathrm{d}(R^2)}{\mathrm{d}S_\mathrm{f}}=0 dSfd(R2)=0(5a)
即有
d(R2)dSf=Lxy2LxxLyy(2LxydLxydSf−1LxxdLxxdSf)(5b)\tag{5b} \frac{\mathrm{d}(R^2)}{\mathrm{d}S_\mathrm{f}}=\frac{L^2_{xy}}{L_{xx}L_{yy}}\bigg(\frac{2}{L_{xy}}\frac{\mathrm{d}L_{xy}}{\mathrm{d}S_\mathrm{f}}-\frac{1}{L_{xx}}\frac{\mathrm{d}L_{xx}}{\mathrm{d}S_\mathrm{f}}\bigg) dSfd(R2)=LxxLyyLxy2(Lxy2dSfdLxy−Lxx1dSfdLxx)(5b)
其中
dLxydSf=−1ln⁡10[∑i=1nyiSi−Sf−1n(∑i=1nyi)(∑i=1n1Si−Sf)]=−1ln⁡10Ly0dLxxdSf=−2ln⁡10[∑i=1nxiSi−Sf−1n(∑i=1nxi)(∑i=1n1Si−Sf)]=−2ln⁡10Lx0(5c)\tag{5c} \frac{\mathrm{d}L_{xy}}{\mathrm{d}S_\mathrm{f}}=\frac{-1}{\ln 10}\bigg[\sum_{i=1}^n\frac{y_i}{S_i-S_\mathrm{f}}-\frac{1}{n}\Big(\sum_{i=1}^ny_i\Big)\Big(\sum_{i=1}^n\frac{1}{S_i-S_\mathrm{f}}\Big)\bigg]=\frac{-1}{\ln 10}L_{y0}\\[8pt] \frac{\mathrm{d}L_{xx}}{\mathrm{d}S_\mathrm{f}}=\frac{-2}{\ln 10}\bigg[\sum_{i=1}^n\frac{x_i}{S_i-S_\mathrm{f}}-\frac{1}{n}\Big(\sum_{i=1}^nx_i\Big)\Big(\sum_{i=1}^n\frac{1}{S_i-S_\mathrm{f}}\Big)\bigg]=\frac{-2}{\ln 10}L_{x0} dSfdLxy=ln10−1[i=1∑nSi−Sfyi−n1(i=1∑nyi)(i=1∑nSi−Sf1)]=ln10−1Ly0dSfdLxx=ln10−2[i=1∑nSi−Sfxi−n1(i=1∑nxi)(i=1∑nSi−Sf1)]=ln10−2Lx0(5c)

将式(5b)、(5c)代入式(5a)，可得
H(Sf)=Ly0Lxy−Lx0Lxx=0(5d)\tag{5d} H(S_\mathrm{f})=\frac{L_{y0}}{L_{xy}}-\frac{L_{x0}}{L_{xx}}=0 H(Sf)=LxyLy0−LxxLx0=0(5d)
于是只需求解非线性方程 H(Sf)H(S_\mathrm{f})H(Sf) 便可得到 SfS_\mathrm{f}Sf ，进一步将 SfS_\mathrm{f}Sf 代入式(4)可得 aaa 、bbb ，则有
{C=10am=−b(6)\tag{6} \begin{cases} C=10^{a}\\[5pt] m=-b \end{cases} ⎩⎨⎧C=10am=−b(6)

计算程序

% 三参数线性拟合
% (S-Sf)^m * N = C
% INPUT
% data_S     : 应力值, 列向量
% data_N     : 疲劳寿命值, 列向量
% k          : 初始点系数(初始点会影响拟合结果), 0~1之间的标量
% OUTPUT
% parameters : [Sf ; m ; C]
function parameters=SN_fit01(data_S,data_N,k)num=length(data_S);data_S=reshape(data_S,num,1);data_N=reshape(data_N,num,1);syms Sff realx=log10(data_S-Sff);x_mean=mean(x);y=log10(data_N);y_mean=mean(y);Lxx=sum(x.^2)-x_mean^2*num;Lxy=sum(x.*y)-x_mean*y_mean*num;Lx0=sum(x./(data_S-Sff))-x_mean*sum(1./(data_S-Sff));Ly0=sum(y./(data_S-Sff))-y_mean*sum(1./(data_S-Sff));H=Lx0/Lxx-Ly0/Lxy;fun_H=matlabFunction(H);fun_m=matlabFunction(-Lxy/Lxx);fun_C=matlabFunction(10^(y_mean-x_mean*Lxy/Lxx));ans_Sf=fzero(fun_H,min(data_S)*k);ans_m=fun_m(ans_Sf);ans_C=fun_C(ans_Sf);parameters=[ans_Sf;ans_m;ans_C];
end

算例

输入：

data_S=[160,120,100,85];
data_N=[96069,273147,434362,2005597];
k=0.8;
y=SN_fit01(data_S,data_N,k);

计算结果：

y=[78.61476407607871.1578247291662316938195.0512843]

非线性拟合方法

将式(1)中的幂函数形式变形为
S=Sf+a∗Nb∗(7a)\tag{7a} S=S_\mathrm{f}+\frac{a^*}{N^{b^*}} S=Sf+Nb∗a∗(7a)
其中
{b∗=1/ma∗=Cb∗(7b)\tag{7b} \begin{cases} b^*=1/m\\[5pt] a^*=C^{b^*} \end{cases} ⎩⎨⎧b∗=1/ma∗=Cb∗(7b)
根据一系列的 SiS_iSi 和 NiN_iNi，对式(7)进行非线性拟合可得 SfS_\mathrm{f}Sf 、mmm 、CCC 。

计算程序

% 三参数非线性拟合
% (S-Sf)^m * N = C
%  S = Sf + a / N^b, 其中 b = 1/m, a = C^b
% INPUT
% data_S     : 应力值, 列向量
% data_N     : 疲劳寿命值, 列向量
% k          : 初始点系数(初始点会影响拟合结果), 0~1之间的标量
% OUTPUT
% parameters : [Sf ; m ; C]
function parameters=SN_fit02(data_S,data_N,k)num=length(data_S);data_S=reshape(data_S,num,1);data_N=reshape(data_N,num,1);fo=fitoptions('Method','NonlinearLeastSquares',...'Lower',[0,0,0],...'Upper',[min(data_S),Inf,Inf],...'StartPoint',[min(data_S)*k,1e4,1]);ft=fittype('Sf+a./N.^b','independent','N',...'dependent','S','options',fo);fa=fit(data_N,data_S,ft);Sf=fa.Sf;m=1/fa.b;C=fa.a^m;parameters=[Sf;m;C];
end

算例

输入：

data_S=[160,120,100,85];
data_N=[96069,273147,434362,2005597];
k=0.8;
y=SN_fit02(data_S,data_N,k);

计算结果：

y=[72.81012886877161.4792077603623571844845.3819234]