三参数 S-N 曲线拟合

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

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

线性拟合方法

将式(1)取对数，可得

\[ \tag{2} m\lg(S-S_\mathrm{f})+\lg N=\lg C \] 令 \(x=\lg(S-S_\mathrm{f})\) 、\(y=\lg N\) 、\(a=\lg C\) 、\(b=-m\) ，可得 \[ \tag{3} y=a+bx \] 对于一系列的应力值 \(S_i\) 、疲劳寿命值 \(N_i\) (\(i=1,2,3,\cdots,n\))，由最小二乘法可得 \[ \tag{4a} \begin{cases} b=L_{xy}/L_{xx}\\[5pt] a=\bar{y}-\bar{x}b \end{cases} \] 线性相关系数 \(R\) 的平方为 \[ \tag{4b} R^2=\frac{L_{xy}^2}{L_{xx}L_{yy}} \] 其中 \[ 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) \] 由于 \(S_\mathrm{f}\) 未知，故式(4)中的表达式均可看作关于 \(S_\mathrm{f}\) 的函数。为了得到最佳的线性相关性，即线性相关系数 \(R\) 的绝对值最大，于是有 \[ \tag{5a} \frac{\mathrm{d}(R^2)}{\mathrm{d}S_\mathrm{f}}=0 \] 即有 \[ \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) \] 其中 \[ \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} \]

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

---

计算程序

% 三参数线性拟合
% (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 real
    x=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.6147640760787
   1.15782472916623
   16938195.0512843]

---

非线性拟合方法

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

---

计算程序

% 三参数非线性拟合
% (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.8101288687716
   1.47920776036235
   71844845.3819234]

---