超空泡模型的力和力矩计算公式

1. 空化器上的阻力和升力

1.1 阻力系数CnentryC_{n_{\text{entry}}}

  • 参数说明：
    -Γ\Gamma：攻角，表示物体与水流方向的夹角。
    -Δ\Delta：空化器倾角，表示空化器的倾斜角度。
    -Ψ0\Psi_0：初始俯仰角，用于描述物体的初始姿态。
    -XnX_n：当前模型的xx坐标，用于跟踪物体位置。
    -RnR_n：空化器的无量纲半径，用于相关计算。
    -CnmaxC_{n_{\text{max}}}：最大法向阻力系数，用于阻力计算。
    -Cx0C_{x0}：空化数为 0 时的阻力系数，用于阻力计算。
    -s1,s2s_1, s_2：中间变量，用于计算阻力系数的过渡值。

  • 公式：

    • 当物体入水时（FlagDive = True）：
      • 若Γ+Δ+Ψ0=−90∘\Gamma + \Delta + \Psi_0 = -90^\circ：
        • 当s1=XnRn≤3.0s_1 = \frac{X_n}{R_n} \leq 3.0时：
          Cnentry=Cnmax+(Cx0−Cnmax)⋅s13.0 C_{n_{\text{entry}}} = C_{n_{\text{max}}} + (C_{x0} - C_{n_{\text{max}}}) \cdot \frac{s_1}{3.0}
        • 否则：
          Cnentry=Cx0 C_{n_{\text{entry}}} = C_{x0}
      • 否则：
        • 当s2=−Xnsin⁡(Γ)2.0Rncos⁡(Γ)≤0.65s_2 = -\frac{X_n \sin(\Gamma)}{2.0 R_n \cos(\Gamma)} \leq 0.65时：
          Cnentry=Cnmax⋅s20.65 C_{n_{\text{entry}}} = C_{n_{\text{max}}} \cdot \frac{s_2}{0.65}
        • 当0.65<s2≤1.40.65 < s_2 \leq 1.4时：
          Cnentry=(Cx0−Cnmax)s2+1.4Cnmax−0.65Cx00.75 C_{n_{\text{entry}}} = \frac{(C_{x0} - C_{n_{\text{max}}}) s_2 + 1.4 C_{n_{\text{max}}} - 0.65 C_{x0}}{0.75}
        • 否则：
          Cnentry=Cx0 C_{n_{\text{entry}}} = C_{x0}
    • 当物体未入水时（FlagDive = False）：
      • 当s1=XnRn≤3.0s_1 = \frac{X_n}{R_n} \leq 3.0时：
        Cnentry=Cnmax+(Cx0−Cnmax)⋅s13.0 C_{n_{\text{entry}}} = C_{n_{\text{max}}} + (C_{x0} - C_{n_{\text{max}}}) \cdot \frac{s_1}{3.0}
      • 否则：
        Cnentry=Cx0 C_{n_{\text{entry}}} = C_{x0}

1.2 法向力系数

  • 参数说明：
    -αJend\alpha_{J_{\text{end}}}：攻角记录参量，用于记录攻角的变化。
    -XcX_c：质心的xx坐标，用于描述物体的质心位置。

  • 公式：
    -xx方向的法向力系数：
    Cnx=Cnentry⋅cos⁡(αJend+Δ180π)⋅cos⁡(Δ180π) C_{nx} = C_{n_{\text{entry}}} \cdot \cos\left(\alpha_{J_{\text{end}}} + \frac{\Delta}{180} \pi\right) \cdot \cos\left(\frac{\Delta}{180} \pi\right)
    -yy方向的法向力系数：
    Cny=−Cnentry⋅cos⁡(αJend+Δ180π)⋅sin⁡(Δ180π) C_{ny} = -C_{n_{\text{entry}}} \cdot \cos\left(\alpha_{J_{\text{end}}} + \frac{\Delta}{180} \pi\right) \cdot \sin\left(\frac{\Delta}{180} \pi\right)

    • 水平力矩：
      Cnm=Xc⋅Cny C_{nm} = X_c \cdot C_{ny}

2. 滑行力和力矩

2.1 下壁面

  • 参数说明：
    -fh1fh1：模型与空泡的下壁面间隙，用于描述空泡特性。
    -FlagCont1FlagCont1：下部第一次接触标识，用于记录接触状态。
    -NcontN_{\text{cont}}：接触次数，用于统计物体接触的次数。
    -CfC_f：表面摩擦系数，用于计算摩擦力。
    -Sw1Sw1：下部浸湿面积，用于描述物体的浸湿状态。
    -B2B^2：空化器的无量纲面积，用于相关计算。
    -Cl2Cl2：上部间隙，用于描述物体与边界的间隙。
    -hshs：修正因子，用于相关计算。
    -f1y,f2yf_{1y}, f_{2y}：中间变量，用于计算滑行力。
    -VpV_p：当前速度，用于描述物体的速度。
    -Yax0Y_{\text{ax0}}：初始空轴线位置，用于描述空泡的初始位置。
    -HXHX：计算步长，用于数值模拟中的递归控制。
    -VyAbsVy_{\text{Abs}}：绝对纵向速度，用于描述物体的纵向速度。
    -Ω\Omega：角速度，用于描述物体的旋转速度。
    -DSDS：空泡变形量/空腔截面的膨胀速度。
    -RcavR_{\text{cav}}：空腔半径数组，用于描述空腔的半径变化。
    -RibTanRibTan：空泡边界速度或模型肋处的切线斜率，用于描述空泡特性。
    -BaseRBaseR：每一节末端的无量纲数，用于描述模型的结构。
    -Imin1I_{\text{min1}}：空泡间隙的位置标识，用于描述空泡间隙的位置。
    -Jc1Jc1：下部接触索引，用于记录接触位置。
    -fdl1fdl1：沾湿面中心或下部浸湿部分长度，用于描述浸湿状态。
    -Arm1Arm1：下方接触力臂，用于计算力矩。

  • 公式：

    • 当下壁面间隙fh1>0.0fh1 > 0.0时：
      • 若未接触过下壁面（FlagCont1 = False），则标记接触并增加接触次数NcontN_{\text{cont}}。
      • 下壁面轴向滑行力系数：
        Csx1=Cf⋅Sw1B2 C_{sx1} = C_f \cdot \frac{Sw1}{B^2}
      • 修正因子hshs：
        hs=−2.0⋅fh1fh1+fh2 hs = -\frac{2.0 \cdot fh1}{fh1 + fh2}
      • 计算f1yf_{1y}和f2yf_{2y}：
        f1y=hs(2.0+hs)(1.0+hs)2 f_{1y} = \frac{hs (2.0 + hs)}{(1.0 + hs)^2}
        f2y=2.0⋅hs1.0+AK⋅hs f_{2y} = \frac{2.0 \cdot hs}{1.0 + AK \cdot hs}
      • 当前速度Vp=VJendV_p = V_{J_{\text{end}}}。
      • 空泡轴的横向速度VycV_{yc}：
        Vyc=Vp⋅Yax0[Jend−Jc1[Imin1]]−Yax0[Jend−Jc1[Imin1]−1]HX V_{yc} = V_p \cdot \frac{Y_{\text{ax0}}[J_{\text{end}} - Jc1[I_{\text{min1}}]] - Y_{\text{ax0}}[J_{\text{end}} - Jc1[I_{\text{min1}}] - 1]}{HX}
      • 模型底部的横向速度V1V_1：
        V1=−VyAbs+Ω⋅(Part[Imin1]−Xc)+Vyc V_1 = -Vy_{\text{Abs}} + \Omega \cdot (Part[I_{\text{min1}}] - X_c) + V_{yc}
      • 空泡边界速度V2V_2：
        V2=−DS[Jc1[Imin1]]2.0⋅π⋅Rcav[Jend,Jc1[Imin1]]+Vp⋅RibTan[Imin1] V_2 = -\frac{DS[Jc1[I_{\text{min1}}]]}{2.0 \cdot \pi \cdot R_{\text{cav}}[J_{\text{end}}, Jc1[I_{\text{min1}}]]} + V_p \cdot RibTan[I_{\text{min1}}]
      • 下壁面横向滑行力系数：
        Csy1=2.0⋅(BaseR[Imin1]Rn)2⋅Vp(V1f1y+V2f2y)VJend2 C_{sy1} = 2.0 \cdot \left(\frac{BaseR[I_{\text{min1}}]}{R_n}\right)^2 \cdot \frac{V_p (V_1 f_{1y} + V_2 f_{2y})}{V_{J_{\text{end}}}^2}
        若Csy1<0.0C_{sy1} < 0.0，则Csy1=0.0C_{sy1} = 0.0。
      • 下方接触力臂：
        Arm1=Part[Imin1]−Xc−0.5⋅fdl1 Arm1 = Part[I_{\text{min1}}] - X_c - 0.5 \cdot fdl1
      • 下壁面接触产生的力矩系数：
        Csm1=−Arm1⋅Csy1 C_{sm1} = -Arm1 \cdot C_{sy1}
    • 当fh1≤0.0fh1 \leq 0.0时：
      • 重置接触标志和相关参数：
        FlagCont1=False,Arm1=0.0,Csx1=0.0,Csy1=0.0,Csm1=0.0 FlagCont1 = False, \quad Arm1 = 0.0, \quad C_{sx1} = 0.0, \quad C_{sy1} = 0.0, \quad C_{sm1} = 0.0

2.2 上壁面

  • 参数说明：
    -fh2fh2：模型与空泡的上壁面间隙，用于描述空泡特性。
    -FlagCont2FlagCont2：上部第一次接触标识，用于记录接触状态。
    -Sw2Sw2：上部浸湿面积，用于描述物体的浸湿状态。
    -Cl1Cl1：下部间隙，用于描述物体与边界的间隙。
    -Imin2I_{\text{min2}}：空泡间隙的位置标识，用于描述空泡间隙的位置。
    -Jc2Jc2：上部接触索引，用于记录接触位置。
    -fdl2fdl2：沾湿面中心或上部浸湿部分长度，用于描述浸湿状态。
    -Arm2Arm2：上方接触力臂，用于计算力矩。

  • 公式：

    • 当上壁面间隙fh2>0.0fh2 > 0.0时：
      • 若未接触过上壁面（FlagCont2 = False），则标记接触并增加接触次数NcontN_{\text{cont}}。
      • 上壁面轴向滑行力系数：
        Csx2=Cf⋅Sw2B2 C_{sx2} = C_f \cdot \frac{Sw2}{B^2}
      • 修正因子hshs：
        hs=−2.0⋅fh2fh1+fh2 hs = -\frac{2.0 \cdot fh2}{fh1 + fh2}
      • 计算f1yf_{1y}和f2yf_{2y}：
        f1y=hs(2.0+hs)(1.0+hs)2 f_{1y} = \frac{hs (2.0 + hs)}{(1.0 + hs)^2}
        f2y=2.0⋅hs1.0+AK⋅hs f_{2y} = \frac{2.0 \cdot hs}{1.0 + AK \cdot hs}
      • 当前速度Vp=VJendV_p = V_{J_{\text{end}}}。
      • 空泡轴的横向速度VycV_{yc}：
        Vyc=Vp⋅Yax0[Jend−Jc2[Imin2]]−Yax0[Jend−Jc2[Imin2]−1]HX V_{yc} = V_p \cdot \frac{Y_{\text{ax0}}[J_{\text{end}} - Jc2[I_{\text{min2}}]] - Y_{\text{ax0}}[J_{\text{end}} - Jc2[I_{\text{min2}}] - 1]}{HX}
      • 模型底部的横向速度V1V_1：
        V1=−VyAbs+Ω⋅(Part[Imin2]−Xc)+Vyc V_1 = -Vy_{\text{Abs}} + \Omega \cdot (Part[I_{\text{min2}}] - X_c) + V_{yc}
      • 空泡边界速度V2V_2：
        V2=DS[Jc2[Imin2]]2.0⋅π⋅Rcav[Jend,Jc2[Imin2]]−Vp⋅RibTan[Imin2] V_2 = \frac{DS[Jc2[I_{\text{min2}}]]}{2.0 \cdot \pi \cdot R_{\text{cav}}[J_{\text{end}}, Jc2[I_{\text{min2}}]]} - V_p \cdot RibTan[I_{\text{min2}}]
      • 上壁面横向滑行力系数：
        Csy2=2.0⋅(BaseR[Imin2]Rn)2⋅Vp(V1f1y+V2f2y)VJend2 C_{sy2} = 2.0 \cdot \left(\frac{BaseR[I_{\text{min2}}]}{R_n}\right)^2 \cdot \frac{V_p (V_1 f_{1y} + V_2 f_{2y})}{V_{J_{\text{end}}}^2}
        若Csy2>0.0C_{sy2} > 0.0，则Csy2=0.0C_{sy2} = 0.0。
      • 上方接触力臂：
        Arm2=Part[Imin2]−Xc−0.5⋅fdl2 Arm2 = Part[I_{\text{min2}}] - X_c - 0.5 \cdot fdl2
      • 上壁面接触产生的力矩系数：
        Csm2=−Arm2⋅Csy2 C_{sm2} = -Arm2 \cdot C_{sy2}
    • 当fh2≤0.0fh2 \leq 0.0时：
      • 重置接触标志和相关参数：
        FlagCont2=False,Arm2=0.0,Csx2=0.0,Csy2=0.0,Csm2=0.0 FlagCont2 = False, \quad Arm2 = 0.0, \quad C_{sx2} = 0.0, \quad C_{sy2} = 0.0, \quad C_{sm2} = 0.0

3. 最终的力和力矩

  • 公式：
    • 总轴向滑行力系数：
      Csx=Csx1+Csx2 C_{sx} = C_{sx1} + C_{sx2}
    • 总横向滑行力系数：
      Csy=Csy1+Csy2 C_{sy} = C_{sy1} + C_{sy2}
    • 总滑行力产生的力矩系数：
      Csm=Csm1+Csm2 C_{sm} = C_{sm1} + C_{sm2}

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
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
def CalcForces(self):
    # CalcForces 子程序是代码中用于计算作用在超空泡模型上的力和力矩的核心部分。
    # 它主要考虑了空化器上的阻力、升力以及模型与空腔壁之间的相互作用力。以下是该子程序的详细分析和解释：
    # CalcForces 的主要功能是计算以下内容：
    # 空化器上的阻力和升力：基于空化数、攻角和模型参数。
    # 滑行力（Planing Forces）：模型与空腔壁接触时产生的力。
    # 力矩：由于滑行力和空化器力产生的力矩。
    # 冲击力：当模型与空腔壁接触时，可能产生的冲击力。
    # 圆周率，用于数学计算
    PI = self.PI
    # 当前计算到的位置，用于跟踪计算进度
    Jend = self.Jend
    # 计算步长，用于数值模拟中的递归控制
    HX = self.HX
    # 空化器的无量纲面积，用于相关计算[^2^]
    B2 = self.B2
    # 空化器的无量纲半径，用于相关计算[^2^]
    Rn = self.Rn
    # 是否入水的标记，用于判断物体是否已进入水中
    FlagDive = self.FlagDive
    # 攻角，表示物体与水流方向的夹角
    Gamma = self.Gamma
    # 空化器倾角，表示空化器的倾斜角度
    Delta = self.Delta
    # 初始俯仰角，用于描述物体的初始姿态
    Psi0 = self.Psi0
    # 控制显示设置，用于决定是否显示某些信息
    SwDisp = self.SwDisp
    # 是否接触的标识，用于判断物体是否与其他物体或边界接触
    SwCont = self.SwCont
    # 计算进度标识，用于跟踪计算的整体进度
    FlagProgress = self.FlagProgress
    # 模型类型，用于区分不同的模型或计算模式
    SwModel = self.SwModel

    # 当前模型的 x 坐标，用于跟踪物体位置
    Xn = self.Xn
    # 最大法向阻力系数，用于阻力计算[^2^]
    CnMax = self.CnMax
    # 空化数为 0 时的阻力系数，用于阻力计算[^2^]
    Cx0 = self.Cx0
    # 攻角记录参量，用于记录攻角的变化[^2^]
    Alpha = self.Alpha
    # 质心的 x 坐标，用于描述物体的质心位置
    Xc = self.Xc
    # 模型与空泡的下壁面间隙，用于描述空泡特性[^2^]
    fh1 = self.fh1
    # 模型与空泡的上壁面间隙，用于描述空泡特性[^2^]
    fh2 = self.fh2
    # 下部第一次接触标识，用于记录接触状态[^2^]
    FlagCont1 = self.FlagCont1
    # 上部第一次接触标识，用于记录接触状态[^2^]
    FlagCont2 = self.FlagCont2
    # 声音标识，用于控制是否发出声音或相关提示[^2^]
    SwSound = self.SwSound
    # 接触次数，用于统计物体接触的次数[^2^]
    Ncont = self.Ncont
    # 表面摩擦系数，用于计算摩擦力[^2^]
    Cf = self.Cf
    # 下部浸湿面积，用于描述物体的浸湿状态[^2^]
    Sw1 = self.Sw1
    # 上部浸湿面积，用于描述物体的浸湿状态[^2^]
    Sw2 = self.Sw2
    # 上部间隙，用于描述物体与边界的间隙[^2^]
    Cl2 = self.Cl2
    # 下部间隙，用于描述物体与边界的间隙[^2^]
    Cl1 = self.Cl1
    # 修正因子，用于相关计算[^2^]
    AK = self.AK
    # 速度记录参量，用于记录速度的变化[^2^]
    V = self.V
    # 初始空轴线位置，用于描述空泡的初始位置[^2^]
    Yax0 = self.Yax0
    # 绝对纵向速度，用于描述物体的纵向速度[^2^]
    VyAbs = self.VyAbs
    # 角速度，用于描述物体的旋转速度[^2^]
    Omega = self.Omega
    # 描述部件位置的数组，用于记录部件的位置信息[^2^]
    Part = self.Part
    # 空泡变形量/空腔截面的膨胀速度
    DS = self.DS
    # 空腔半径数组，用于描述空腔的半径变化[^2^]
    Rcav = self.Rcav
    # 空泡边界速度或模型肋处的切线斜率，用于描述空泡特性[^2^]
    RibTan = self.RibTan
    # 每一节末端的无量纲数，用于描述模型的结构[^2^]
    BaseR = self.BaseR
    # 空泡间隙的位置标识，用于描述空泡间隙的位置[^2^]
    Imin1 = self.Imin1
    # 空泡间隙的位置表示
    Imin2 = self.Imin2
    # 下部接触索引，用于记录接触位置[^2^]
    Jc1 = self.Jc1
    # 上部接触索引，用于记录接触位置[^2^]
    Jc2 = self.Jc2
    # 沾湿面中心或下部浸湿部分长度，用于描述浸湿状态[^2^]
    fdl1 = self.fdl1
    # 沾湿面中心或上部浸湿部分长度，用于描述浸湿状态[^2^]
    fdl2 = self.fdl2

    # 作用在空泡器上的力和力矩的分量
    if FlagDive:
        if Gamma + Delta + Psi0 == -90:
            s1 = Xn / Rn
            if s1 <= 3.0:
                CnEntry = CnMax + (Cx0 - CnMax) * s1 / 3.0
            else:
                CnEntry = Cx0
        else:
            s2 = -Xn * np.sin(Gamma) / (np.cos(Gamma) * 2.0 * Rn)
            if s2 <= 0.65:
                CnEntry = CnMax * s2 / 0.65
            elif s2 <= 1.4:
                CnEntry = ((Cx0 - CnMax) * s2 + 1.4 * CnMax - 0.65 * Cx0) / 0.75
            else:
                CnEntry = Cx0
    else:
        # （这是我们所要关注的重点部分）
        s1 = Xn / Rn
        if s1 <= 3.0:
            CnEntry = CnMax + (Cx0 - CnMax) * s1 / 3.0
        else:
            CnEntry = Cx0
    # x方向的法向力系数
    Cnx = CnEntry * np.cos(Alpha[Jend] + Delta / 180*np.pi()) * np.cos(Delta/180*np.pi())
    # y方向的法向力系数
    Cny = -CnEntry * np.cos(Alpha[Jend] + Delta / 180*np.pi()) * np.sin(Delta/180*np.pi())
    # 水平力矩
    Cnm = Xc * Cny

    # Planing forces and moments calculation
    # 平面力和力矩的分量
    if SwModel == 0:
        if fh1 > 0.0:
            if not FlagCont1:
                FlagCont1 = True
                Ncont += 1
                if not FlagProgress and (SwCont == 1 or SwDisp == 1):
                    print(3)
            # Cf是表面摩擦阻力系数，Sw1是沾湿面积，但具体计算公式推导过程及意义暂不明确。现已探明沾湿面积计算公式
            # 下壁面轴向滑行力系数，用于计算相关受力[^2^]
            Csx1 = Cf * Sw1 / B2
            # fh2是间隙
            fh2 = Cl2[Imin1]  # Imin1 assumed to be 0 for example
            # 这里和第六讲的h一把是等价表述
            # 修正因子hs(CalcForces)
            hs = -2.0 * fh1 / (fh1 + fh2)
            # 见第六讲6.15
            f1y = hs * (2.0 + hs) / (1.0 + hs) ** 2
            # 见第六讲6.15
            f2y = 2.0 * hs / (1.0 + AK * hs)
            # 当前速度
            Vp = V[Jend]
            # 空泡轴的横向速度
            Vyc = Vp * (Yax0[Jend - Jc1[Imin1]] - Yax0[Jend - Jc1[Imin1] - 1]) / HX
            # 见第六讲，是模型底部的横向速度
            V1 = -VyAbs + Omega * (Part[Imin1] - Xc) + Vyc
            # !揣测为空泡边界速度，但表达形式与文献一致或等价性还需要确认
            # !①这里加条件控制语句，来控制是否有冲击力, 但鉴于没有角速度量可以确保在这里被准确无望的使用，或考虑在der中进行定义。
            V2 = -DS[Jc1[Imin1]] / (2.0 * PI * Rcav[Jend, Jc1[Imin1]]) + Vp * RibTan[Imin1]

            # 下壁面横向滑行力系数，用于计算相关受力[^2^]
            Csy1 = 2.0 * (BaseR[Imin1] / Rn) ** 2 * Vp * (V1 * f1y + V2 * f2y) / V[Jend] ** 2
            # 吸引力是不被允许的
            if Csy1 < 0.0:
                Csy1 = 0.0
            # 下方接触力臂 6.16，最后减去的那项是展示面中心的意思
            Arm1 = Part[Imin1] - Xc - 0.5 * fdl1
            # 下壁面接触产生的力矩系数，用于计算相关力矩[^2^]
            Csm1 = -Arm1 * Csy1
        else:
            # 如果要设置常力，可以在这里设置。力的位置只需要改变力臂大小即可，即对CSM1做变动
            FlagCont1 = False
            Arm1 = 0.0
            Csx1 = 0.0
            Csy1 = 0.0
            Csm1 = 0.0

        ################上壁面
        if fh2 > 0.0:
            if not FlagCont2:

                FlagCont2 = True
                Ncont += 1
                if not FlagProgress and (SwCont == 1 or SwDisp == 1):
                    print(3)
            # 上壁面轴向滑行力系数(CalcForces)
            Csx2 = Cf * Sw2 / B2
            fh1 = Cl1[Imin2]  # Imin2 assumed to be 0 for example
            hs = -2.0 * fh2 / (fh1 + fh2)
            f1y = hs * (2.0 + hs) / (1.0 + hs) ** 2
            f2y = 2.0 * hs / (1.0 + AK * hs)
            Vp = V[Jend]
            Vyc = Vp * (Yax0[Jend - Jc2[Imin2]] - Yax0[Jend - Jc2[Imin2] - 1]) / HX
            V1 = -VyAbs + Omega * (Part[Imin2] - Xc) + Vyc
            V2 = DS[Jc2[Imin2]] / (2.0 * PI * Rcav[Jend, Jc2[Imin2]]) - Vp * RibTan[Imin2]
            # 上壁面横向滑行力系数(CalcForces)
            Csy2 = 2.0 * (BaseR[Imin2] / Rn) ** 2 * Vp * (V1 * f1y + V2 * f2y) / V[Jend] ** 2

            if Csy2 > 0.0:
                # 通过这里控制了上表面接触产生的滑行力只能往下，就是说，Csy2如果存在必定小于0。综合可知。力矩顺时针为正，逆时针为负。
                Csy2 = 0.0
            # 上接触力臂
            Arm2 = Part[Imin2] - Xc - 0.5 * fdl2
            # 上表面滑行力产生的力矩是顺时针，是正的。
            # 上壁面接触产生的力矩系数，用于计算相关力矩[^2^]
            Csm2 = -Arm2 * Csy2
        else:
            FlagCont2 = False
            Arm2 = 0.0
            Csx2 = 0.0
            Csy2 = 0.0
            Csm2 = 0.0

    # Final forces and moments
    # 这里好像只要SwModel为0了就会有报错，注意一下
    Csx = Csx1 + Csx2
    Csy = Csy1 + Csy2
    Csm = Csm1 + Csm2

    self.Cnx = Cnx
    self.Cny = Cny
    self.Cnm = Cnm
    self.Csx = Csx
    self.Csy = Csy
    self.Csm = Csm

    self.fh2 = fh2
    self.fh1 = fh1
    self.FlagCont1 = FlagCont1
    self.Ncont = Ncont
    self.FlagCont2 = FlagCont2