CalcStress 函数的数学公式表示

函数描述

CalcStress 函数用于从当前速度计算模型截面的弹性应力。该函数涉及多个物理量和计算步骤，包括轴向内力、剪切内力、正应力和剪切应力的计算。

参数说明

  • VcurV_{\text{cur}}：当前速度
  • V0V_0：初始速度
  • SnS_n：截面面积
  • Cn,X2strC_{n,\text{X2str}}、Cs,X2strC_{s,\text{X2str}}、SG,X2strS_{G,\text{X2str}}：应力相关参数
  • SGS_G：空化数数组
  • CnxC_{nx}、CsyC_{sy}：力系数
  • FlagStress\text{FlagStress}：应力计算标志
  • JendJ_{\text{end}}：应力计算结束索引
  • XmstrX_{\text{mstr}}：截面位置
  • LmmL_{\text{mm}}：模型长度（毫米）
  • LfL_f：前体长度
  • ρf\rho_f、ρa\rho_a：流体和空气密度
  • SmS_m：截面面积
  • TMkgTM_{\text{kg}}：模型总质量
  • XcX_c：质心位置
  • IcI_c：惯性矩
  • RmstrR_{\text{mstr}}：截面半径
  • RnmmR_{\text{nmm}}：归一化半径
  • CoBeg1\text{CoBeg1}、CoEnd1\text{CoEnd1}、CoFact1\text{CoFact1}：应力集中因子（前体）
  • CoBeg2\text{CoBeg2}、CoEnd2\text{CoEnd2}、CoFact2\text{CoFact2}：应力集中因子（后体）
  • FlagStress\text{FlagStress}：截面应力计算标志
  • X1,strX_{1,\text{str}}：特定截面位置

数学公式表示

1. 空化数和力系数的选择

根据 FlagStress 的值选择空化数和力系数：
{SG,cur=SG,X2strCn,cur=Cn,X2strCs,cur=∣Cs,X2str∣if FlagStress=TrueSG,cur=SG[Jend]Cn,cur=CnxCs,cur=∣Csy∣otherwise \begin{cases} S_{G,\text{cur}} = S_{G,\text{X2str}} \\ C_{n,\text{cur}} = C_{n,\text{X2str}} \\ C_{s,\text{cur}} = |C_{s,\text{X2str}}| & \text{if } \text{FlagStress} = \text{True} \\ S_{G,\text{cur}} = S_G[J_{\text{end}}] \\ C_{n,\text{cur}} = C_{nx} \\ C_{s,\text{cur}} = |C_{sy}| & \text{otherwise} \end{cases}

2. 轴向内力NxN_x的计算

Nx[0]=103⋅(Vcur⋅V0)2⋅Sn⋅Cn,cur⋅(1+SG,cur)/2(单位：牛顿) N_x[0] = 10^3 \cdot (V_{\text{cur}} \cdot V_0)^2 \cdot S_n \cdot C_{n,\text{cur}} \cdot \left(1 + S_{G,\text{cur}}\right) / 2 \quad \text{(单位：牛顿)}
Ax=Nx[0]TMkg A_x = \frac{N_x[0]}{TM_{\text{kg}}}
Nx[i]=Nx[i−1]−103⋅Ax⋅(0.5⋅0.005⋅Lm⋅(ρ1⋅Sm[i−1]+ρ2⋅Sm[i]))for i=1,2,…,200 N_x[i] = N_x[i-1] - 10^3 \cdot A_x \cdot \left(0.5 \cdot 0.005 \cdot L_m \cdot (\rho_1 \cdot S_m[i-1] + \rho_2 \cdot S_m[i])\right) \quad \text{for } i = 1, 2, \ldots, 200
其中：
ρ1={ρfif Xmstr[i−1]⋅Lmm<Lfρaotherwise \rho_1 = \begin{cases} \rho_f & \text{if } X_{\text{mstr}}[i-1] \cdot L_{\text{mm}} < L_f \\ \rho_a & \text{otherwise} \end{cases}
ρ2={ρfif Xmstr[i]⋅Lmm<Lfρaotherwise \rho_2 = \begin{cases} \rho_f & \text{if } X_{\text{mstr}}[i] \cdot L_{\text{mm}} < L_f \\ \rho_a & \text{otherwise} \end{cases}

3. 正应力σx1\sigma_{x1}的计算

σx1=10−6⋅NxSm(单位：MPa) \sigma_{x1} = 10^{-6} \cdot \frac{N_x}{S_m} \quad \text{(单位：MPa)}

4. 剪切内力QyQ_y的计算

Fs=103⋅(Vcur⋅V0)2⋅Sn⋅Cs,cur/2(单位：牛顿) F_s = 10^3 \cdot (V_{\text{cur}} \cdot V_0)^2 \cdot S_n \cdot C_{s,\text{cur}} / 2 \quad \text{(单位：牛顿)}
Ay=FsTMkg A_y = \frac{F_s}{TM_{\text{kg}}}
Qy[i]=103⋅Ay⋅(0.5⋅0.005⋅Lm⋅(ρ1⋅Sm[i−1]+ρ2⋅Sm[i]))for i=1,2,…,200 Q_y[i] = 10^3 \cdot A_y \cdot \left(0.5 \cdot 0.005 \cdot L_m \cdot (\rho_1 \cdot S_m[i-1] + \rho_2 \cdot S_m[i])\right) \quad \text{for } i = 1, 2, \ldots, 200

5. 剪切应力τy\tau_y的计算

τy=10−6⋅4⋅Qy3⋅Sm(单位：MPa) \tau_y = 10^{-6} \cdot \frac{4 \cdot Q_y}{3 \cdot S_m} \quad \text{(单位：MPa)}

6. 弯矩MzM_z的计算

Mz[i]=Xmstr[i]⋅Lm⋅Qy[i]−(Ay−Fs⋅Lm⋅(1−Xc)Ic)⋅103⋅(0.5⋅0.005⋅Lm2⋅(ρ1⋅Xmstr[i−1]⋅Sm[i−1]+ρ2⋅Xmstr[i]⋅Sm[i]))for i=1,2,…,200 M_z[i] = X_{\text{mstr}}[i] \cdot L_m \cdot Q_y[i] - \left(A_y - \frac{F_s \cdot L_m \cdot (1 - X_c)}{I_c}\right) \cdot 10^3 \cdot \left(0.5 \cdot 0.005 \cdot L_m^2 \cdot (\rho_1 \cdot X_{\text{mstr}}[i-1] \cdot S_m[i-1] + \rho_2 \cdot X_{\text{mstr}}[i] \cdot S_m[i])\right) \quad \text{for } i = 1, 2, \ldots, 200

7. 正应力σx2\sigma_{x2}的计算

σx2[i]=10−6⋅Mz[i]⋅Rmstr[i]⋅Rnmm/103π⋅(Rmstr[i]⋅Rnmm/103)4/4(单位：MPa) \sigma_{x2}[i] = 10^{-6} \cdot \frac{M_z[i] \cdot R_{\text{mstr}}[i] \cdot R_{\text{nmm}} / 10^3}{\pi \cdot (R_{\text{mstr}}[i] \cdot R_{\text{nmm}} / 10^3)^4 / 4} \quad \text{(单位：MPa)}

8. 总正应力σx\sigma_x的计算

σx=σx1+σx2(单位：MPa) \sigma_x = \sigma_{x1} + \sigma_{x2} \quad \text{(单位：MPa)}

9. 应力集中处理

σx[i]={σx[i]⋅CoFact1if CoBeg1<Xmstr[i]⋅Lmm<CoEnd1σx[i]⋅CoFact2if CoBeg2<Xmstr[i]⋅Lmm<CoEnd2σx[i]otherwise \sigma_x[i] = \begin{cases} \sigma_x[i] \cdot \text{CoFact1} & \text{if } \text{CoBeg1} < X_{\text{mstr}}[i] \cdot L_{\text{mm}} < \text{CoEnd1} \\ \sigma_x[i] \cdot \text{CoFact2} & \text{if } \text{CoBeg2} < X_{\text{mstr}}[i] \cdot L_{\text{mm}} < \text{CoEnd2} \\ \sigma_x[i] & \text{otherwise} \end{cases}

10. 最大正应力和最大剪切应力的计算

σmax=max⁡(σx) \sigma_{\text{max}} = \max(\sigma_x)
Xσmax=Xmstr[argmax(σx)] X_{\sigma_{\text{max}}} = X_{\text{mstr}}[\text{argmax}(\sigma_x)]
τmax=max⁡(∣τy∣) \tau_{\text{max}} = \max(|\tau_y|)
Xτmax=Xmstr[argmax(∣τy∣)] X_{\tau_{\text{max}}} = X_{\text{mstr}}[\text{argmax}(|\tau_y|)]

11. 特定截面的应力值

如果需要，计算特定截面的应力值：
QX1,str=interp(X1,str,Xmstr,Qy) Q_{X_{1,\text{str}}} = \text{interp}(X_{1,\text{str}}, X_{\text{mstr}}, Q_y)
τX1,str=interp(X1,str,Xmstr,τy) \tau_{X_{1,\text{str}}} = \text{interp}(X_{1,\text{str}}, X_{\text{mstr}}, \tau_y)
MX1,str=interp(X1,str,Xmstr,Mz) M_{X_{1,\text{str}}} = \text{interp}(X_{1,\text{str}}, X_{\text{mstr}}, M_z)
σX1,str=interp(X1,str,Xmstr,σx) \sigma_{X_{1,\text{str}}} = \text{interp}(X_{1,\text{str}}, X_{\text{mstr}}, \sigma_x)

Python代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
def CalcStress(self, Vcur):
    # , Vcur, V0, Sn, Cn_X2str, Cs_X2str, SG_X2str, SG, Cnx, Csy, FlagStress, Jend,
    #                     Xmstr, Lmm, Lf, Rhof, Rhoa, Sm, TMkg, Xc, Ic, Rmstr, Rnmm, CoBeg1, CoEnd1, CoFact1,
    #                     CoBeg2, CoEnd2, CoFact2, FlagStress_section, X1str
    """
    从当前速度计算模型截面的弹性应力。
    参数：
    - Vcur: 当前速度
    - V0: 初始速度
    - Sn: 截面面积
    - Cn_X2str, Cs_X2str, SG_X2str: 应力相关参数
    - SG: 空化数数组
    - Cnx, Csy: 力系数
    - FlagStress: 应力计算标志
    - Jend: 应力计算结束索引
    - Xmstr: 截面位置
    - Lmm: 模型长度
    - Lf: 前体长度
    - Rhof, Rhoa: 流体和空气密度
    - Sm: 截面面积
    - TMkg: 模型总质量
    - Xc: 质心位置
    - Ic: 惯性矩
    - Rmstr: 截面半径
    - Rnmm: 归一化半径
    - CoBeg1, CoEnd1, CoFact1: 应力集中因子
    - CoBeg2, CoEnd2, CoFact2: 应力集中因子
    - FlagStress_section: 截面应力计算标志
    - X1str: 特定截面位置
    """
    # 应力计算标志，这一版的好像就是false
    FlagStress = self.FlagStress
    # 运动应力菜单对应的，这里好像是交叉引用了
    # 好像是从SG过来的，
    # SGstr(Istr) = SG(Jend)
    # SG_X2str = SGstr(1)（这里是插值）
    SG_X2str = self.SG_X2str
    Cn_X2str = self.Cn_X2str
    Cs_X2str = self.Cs_X2str
    # 空化数存储合集？（会在Dynamic中定义）
    SG = self.SG
    # 当前计算位置
    Jend = self.Jend
    # 几个力矩系数
    Cnx = self.Cnx
    Csy = self.Csy
    # 初始速度
    V0 = self.V0
    # 空化器面积
    Sn = self.Sn
    # 总质量
    TMkg = self.TMkg
    # 模型前体（Forebody）部分的材料密度（kg/m^3）
    Rhof = self.Rhof
    # 模型半径分布的x坐标，用于描述模型的半径变化
    Xmstr = self.Xmstr
    # 模型总长度，单位毫米(CalcParameters)
    Lmm = self.Lmm
    # 模型前体（Forebody）的长度
    Lf = self.Lf
    # 模型长度，单位米(CalcParameters)
    Lm = self.Lm
    # Rhoa (float): 后体材料密度（kg/m^3）。
    Rhoa = self.Rhoa
    # 模型半径分布的横截面积，用于描述模型的几何特性(CalcParameters)
    Sm = self.Sm
    # 质心的X坐标(CalcModel定义)
    Xc = self.Xc
    # 模型质心转动惯量(CalcModel)
    Ic = self.Ic
    # 存储了模型在不同位置的无量纲半径值(ModelRadii)
    Rmstr = self.Rmstr
    # 空化器半径，单位毫米(CalcParameters)
    Rnmm = self.Rnmm
    # # 应力集中定义参量
    # 前体开始间隔
    CoBeg1 = self.CoBeg1
    # 前体结束间隔
    CoEnd1 = self.CoEnd1
    # 前体增益系数
    CoFact1 = self.CoFact1
    # 后体开始间隔，单位毫米
    CoBeg2 = self.CoBeg2
    # 后体结束间隔
    CoEnd2 = self.CoEnd2
    # 后体增益系数
    CoFact2 = self.CoFact2
    # 应力计算参量，想放弃
    X1str = self.X1str

    if FlagStress:
        # 如果是“运动/应力”菜单功能，则使用特定的应力参数
        SGcur = SG_X2str
        Cncur = Cn_X2str
        Cscur = abs(Cs_X2str)
    else:
        # 否则使用动态计算中的当前值
        SGcur = SG[Jend]
        Cncur = Cnx
        Cscur = abs(Csy)

    # 计算轴向内力（与SCAV相同）
    Nx = np.zeros(201)
    Nx[0] = 1e3 * (Vcur * V0) ** 2 * Sn * Cncur * (1 + SGcur) / 2  # 单位：牛顿
    Ax = Nx[0] / TMkg  # 模型纵向加速度，单位：m/s^2

    for i in range(1, 201):
        Rho1 = Rhof if Xmstr[i - 1] * Lmm < Lf else Rhoa
        Rho2 = Rhof if Xmstr[i] * Lmm < Lf else Rhoa
        work = 0.5 * 0.005 * Lm * (Rho1 * Sm[i - 1] + Rho2 * Sm[i])  # 梯形积分法
        Nx[i] = Nx[i - 1] - 1e3 * Ax * work

    # 计算由轴向力Fn引起的正应力
    SGx1 = 1e-6 * Nx / Sm  # 单位：MPa

    # 计算剪切内力Qy
    Fs = 1e3 * (Vcur * V0) ** 2 * Sn * Cscur / 2  # 单位：牛顿

    Qy = np.zeros(201)
    Ay = Fs / TMkg  # 模型横向加速度，单位：m/s^2
    Int1 = 0.0

    for i in range(1, 201):
        Rho1 = Rhof if Xmstr[i - 1] * Lmm < Lf else Rhoa
        Rho2 = Rhof if Xmstr[i] * Lmm < Lf else Rhoa
        Int1 += 0.5 * 0.005 * Lm * (Rho1 * Sm[i - 1] + Rho2 * Sm[i])  # 梯形积分法
        Qy[i] = 1e3 * Ay * Int1

    # 计算剪切应力
    TAUy = 1e-6 * 4 * Qy / (3 * Sm)  # 单位：MPa

    # 计算弯矩Mz
    Mz = np.zeros(201)
    Int2 = 0.0
    Int3 = 0.0
    work1 = Fs * Lm * (1 - Xc) / Ic

    for i in range(1, 201):
        Rho1 = Rhof if Xmstr[i - 1] * Lmm < Lf else Rhoa
        Rho2 = Rhof if Xmstr[i] * Lmm < Lf else Rhoa
        Int2 += 0.5 * 0.005 * Lm ** 2 * (Rho1 * Xmstr[i - 1] * Sm[i - 1] + Rho2 * Xmstr[i] * Sm[i])
        Int3 += 0.5 * 0.005 * Lm ** 3 * (Rho1 * Xmstr[i - 1] ** 2 * Sm[i - 1] + Rho2 * Xmstr[i] ** 2 * Sm[i])
        Mz[i] = Xmstr[i] * Lm * Qy[i] - (Ay - work1 * Xc * Lm) * 1e3 * Int2 - work1 * 1e3 * Int3

    # 计算由横向力Fs引起的正应力
    SGx2 = np.zeros(201)
    for i in range(201):
        Rmm = Rmstr[i] * Rnmm / 1e3  # 单位：米
        Jz = np.pi * Rmm ** 4 / 4
        SGx2[i] = 1e-6 * Mz[i] * Rmm / Jz  # 单位：MPa

    # 计算总正应力
    SGx = SGx1 + SGx2  # 单位：MPa

    # 应力集中处理
    for i in range(201):
        if CoBeg1 < Xmstr[i] * Lmm < CoEnd1:
            SGx[i] *= CoFact1
        elif CoBeg2 < Xmstr[i] * Lmm < CoEnd2:
            SGx[i] *= CoFact2

    # 计算最大正应力
    SGmax = np.max(SGx)
    X_SGmax = Xmstr[np.argmax(SGx)]

    # 计算最大剪切应力
    TAUmax = np.max(np.abs(TAUy))
    X_TAUmax = Xmstr[np.argmax(np.abs(TAUy))]

    # 如果需要，计算特定截面的应力值
    if FlagStress:
        Q_X1str = np.interp(X1str, Xmstr, Qy)
        TAU_X1str = np.interp(X1str, Xmstr, TAUy)
        M_X1str = np.interp(X1str, Xmstr, Mz)
        SG_X1str = np.interp(X1str, Xmstr, SGx)

    # 返回参量
    # 轴向内力（与SCAV相同）
    self.Nx = Nx
    # 计算由轴向力Fn引起的正应力
    self.SGx1 = SGx1
    # 剪切内力？横向剪切力？
    self.Qy = Qy
    # 剪切应力
    self.TAUy = TAUy
    # 弯矩
    self.Mz = Mz
    # 计算由横向力Fs引起的正应力
    self.SGx2 = SGx2
    # 总应力
    self.SGx = SGx
    # 最大正应力
    self.SGmax = SGmax
    # 最大正应力位置
    self.X_SGmax = X_SGmax
    # 最大切应力
    self.TAUmax = TAUmax
    # 最大切应力位置
    self.X_TAUmax = X_TAUmax
    # 用户指定界面剪切力
    self.Q_X1str = Q_X1str
    # 用户指定界面剪切应力
    self.TAU_X1str = TAU_X1str
    # 指定截面的弯矩
    self.M_X1str = M_X1str
    # 指定界面的正应力
    self.SG_X1str = SG_X1str