奇异函数解梁的剪力和弯矩

零件工程师

用Mathcad Prime3定义了个奇异函数，相较于Excel，一旦定义了函数，针对不同项目的具体支反力、力矩，查表后快速写入弯矩、剪力函数。比Excel要方便些。

Fluid

零件工程师

Fluid 发表于 2023-10-9 16:07

大侠那个c矩阵转置，看得出来是老手，不过Mathcad这软件想把结果自动生成一列还真没啥好方法，只能这样转置。

Fluid

零件工程师 发表于 2023-10-9 16:43
大侠那个c矩阵转置，看得出来是老手，不过Mathcad这软件想把结果自动生成一列还真没啥好方法，只能这样转 ...

大侠我之前说过mathcad是你教我的, 我计算能力很差, 只能用用软件了.

可以这样, 就直接将结果指定给未知数了. 但必须是结果出来了,在等式的最前面加 := , 再把未知数矩阵写在:=前面,等式成立.

零件工程师

Fluid 发表于 2023-10-9 16:50
大侠我之前说过mathcad是你教我的, 我计算能力很差, 只能用用软件了.

可以这样, 就直接将结果指定给未知 ...

大侠，我用Python写了个梁的计算，可以试试，加载参数没有变，输出的剪力图没问题，弯矩图在突变处没有体现，留个小尾巴，估计是函数定义的时候没做过多区分，有时间再搞。

1. 
   
2. # 引入外部库，创建程序环境
   
3. import numpy as np  # 引入numpy并简称为np
   
4. import matplotlib.pyplot as plt  # 引入matplotlib环境并简称为plt
   
5. 
   
6. # 定义剪力方程（基于奇异函数规则）
   
7. # 定义剪力方程，引入变量 x：梁的横坐标；a:外部载荷位置：n:函数阶次
   
8. def shear_force(x, a, n):  
   
9.     if n < 0:
   
10.         if x == a:
    
11.             return 1
    
12.         elif x != a:
    
13.             return 0
    
14.     elif n >= 0:
    
15.         if x < a:
    
16.             return 0
    
17.         elif x >= a:
    
18.             return P_Shear * (x - a) ** n
    
19. 
    
20. # 定义弯矩方程（基于奇异函数规则）
    
21. # 定义剪弯矩程，引入变量 x：梁的横坐标；a:外部载荷位置：n:函数阶次
    
22. def moment_value(x, a, n):
    
23.     if n < 0:
    
24.         if x == a:
    
25.             return 1
    
26.         elif x != a:
    
27.             return 0
    
28.     elif n >= 0:
    
29.         if x < a:
    
30.             return 0
    
31.         elif x >= a:
    
32.             return P_Moment * (x - a) ** ( n + 1 )        
    
33.         
    
34. # 引入用于奇异函数计算时所需的系数矩阵
    
35. 
    
36. # 创建外部载荷矩阵，数量、数值可根据具体案例修改
    
37. # 用于剪力计算的载荷系数
    
38. matrix_Load_For_ShearForce = np.array([73.75, -120, -45, -30, 30, 61.25])
    
39. # 用于弯矩计算的载荷系数，注意根据奇异函数规则表更新系数，与剪力的载荷系数有所不同
    
40. matrix_Load_For_Moment = np.array([73.75, -120, -45, -15, 15, 61.25])
    
41. 
    
42. # 创建外部载荷位置矩阵，数量、数值可根据具体案例修改
    
43. matrix_Load_Loaction = np.array([0, 2, 4, 6, 9, 12])
    
44. 
    
45. # 创建阶次矩阵，数量、数值可根据具体案例修改
    
46. matrix_n = np.array([0, -1, 0, 1, 1, 0])
    
47. 
    
48. # 定义横坐标x的起点、终点、细分数量,细分数量影响计算精度
    
49. x = np.linspace(0, 12, 5000)
    
50. # 为shear_force_values生成与x矩阵相同的全零数组
    
51. shear_force_values = np.zeros_like(x)
    
52. 
    
53. # 为moment_values生成与x矩阵相同的全零数组
    
54. moment_values = np.zeros_like(x)
    
55. 
    
56. # 开始遍历矩阵
    
57. for i in range(len(matrix_Load_Loaction)):
    
58.     P_Shear = matrix_Load_For_ShearForce[i] # 创建P并将负载矩阵中的序号为i元素依次赋值给P
    
59.     P_Moment = matrix_Load_For_Moment[i] # 创建P并将负载矩阵中的序号为i元素依次赋值给P
    
60.     a = matrix_Load_Loaction[i] # 创建a并将负载位置矩阵中的序号为i元素依次赋值给a
    
61.     n = matrix_n[i] # 创建n并将阶次矩阵中的序号为i元素依次赋值给n
    
62. 
    
63. # 将shear_force的函数变成一个可以接受数组作为输入的向量化函数shear_force_values  
    
64. # 将x、a和n作为参数传递给向量化的shear_force函数。
    
65. # 向量化函数将对每个元素进行计算，并返回一个具有相同形状的数组，其中包含计算后的结果
    
66. # shear_force_values +=将上述乘积添加到名为shear_force_values的变量中
    
67. # 如果shear_force_values之前已经有值，这将进行逐元素的相加。
    
68. # 这样，shear_force_values逐渐累积了每个循环迭代中的剪力值
    
69.     shear_force_values += np.vectorize(shear_force)(x, a, n)
    
70.    
    
71. # 同上
    
72.     moment_values += np.vectorize(moment_value)(x, a, n)
    
73. 
    
74. # 这行代码的目的是计算剪力值，并将结果累积在变量shear_force_values中
    
75. # 通过使用向量化函数和广播机制，可以高效地对输入数组的每个元素进行计算
    
76. # 并利用向量化的特性，避免了显式地编写循环来处理数组的每个元素
    
77.    
    
78. def plot_shear_force_diagram(x, shear_force_values):
    
79.     plt.plot(x, shear_force_values)
    
80.     plt.xlabel('Position (m)')
    
81.     plt.ylabel('Shear Sorce (kN)')
    
82.     plt.title('Shear Force Diagram')
    
83.     plt.grid(True)
    
84. 
    
85. def plot_moment_diagram(x, moment_values):
    
86.     plt.plot(x, moment_values)
    
87.     plt.xlabel('Position (m)')
    
88.     plt.ylabel('Moment Values (kN·m)')
    
89.     plt.title('Moment Diagram')
    
90.     plt.grid(True)
    
91.    
    
92. # 创建一个包含两个子图的图形
    
93. fig, axs = plt.subplots(2)
    
94. 
    
95. # 在第一个子图中绘制剪力图
    
96. plt.sca(axs[0])
    
97. plot_shear_force_diagram(x, shear_force_values)
    
98. plt.ylim(-100, 100)
    
99. plt.yticks(np.arange(0, 120, 40))
    
100. 
     
101. # 在第二个子图中绘制弯矩图
     
102. plt.sca(axs[1])
     
103. plot_moment_diagram(x, moment_values)
     
104. plt.ylim(-100, 100)
     
105. plt.yticks(np.arange(0, 500, 100))
     
106. 
     
107. # 调整子图之间的间距
     
108. plt.subplots_adjust(hspace=1)
     
109. 
     
110. # 显示图形
     
111. plt.grid(True)
     
112. plt.show()

复制代码

零件工程师

1. import math
   
2. import sympy
   
3. 
   
4. # 将角度转换为弧度
   
5. angle_a = 55  # 输入角度a的角度制角度
   
6. angle_b = 40  # 输入角度b的角度制角度
   
7. r_a = math.radians(angle_a)  # 角度制转换为弧度制
   
8. r_b = math.radians(angle_b)  # 角度制转换为弧度制
   
9. 
   
10. # 列平衡方程
    
11. F_1, F_2, P, = sympy.symbols("F_1, F_2, P")  # 定义符号计算用的变量
    
12. eq1 = F_2 * math.cos(r_a) - F_1 * math.sin(r_b)  # 平衡方程1
    
13. eq2 = F_2 * math.sin(r_a) + F_1 * math.cos(r_b) - P # 平衡方程2
    
14. s = sympy.solve([eq1, eq2], [F_1,F_2], dict=True)  # 求解平衡方程组
    
15. print(s)  #打印结果

复制代码

零件工程师

Fluid 发表于 2023-10-9 16:50
大侠我之前说过mathcad是你教我的, 我计算能力很差, 只能用用软件了.

可以这样, 就直接将结果指定给未知 ...

问题解决了，弯矩方程的奇异函数规则设置的有问题，应该评判n + 1是否小于0，因为剪力方程的规则中已经用过一次n了，会返回值。

1. # 引入外部库，创建程序环境
   
2. import numpy as np  # 引入numpy并简称为np
   
3. import matplotlib.pyplot as plt  # 引入matplotlib环境并简称为plt
   
4. 
   
5. # 定义剪力方程（基于奇异函数规则）
   
6. # 定义剪力方程，引入变量 x：梁的横坐标；a:外部载荷位置：n:函数阶次
   
7. def shear_force(x, a, n):  
   
8.     if n < 0:
   
9.         if x == a:
   
10.             return 1
    
11.         elif x != a:
    
12.             return 0
    
13.     elif n >= 0:
    
14.         if x < a:
    
15.             return 0
    
16.         elif x >= a:
    
17.             return P_Shear * (x - a) ** n
    
18. 
    
19. # 定义弯矩方程（基于奇异函数规则）
    
20. # 定义剪弯矩程，引入变量 x：梁的横坐标；a:外部载荷位置：n:函数阶次
    
21. def moment_value(x, a, n):
    
22.     if n + 1 < 0:
    
23.         if x == a:
    
24.             return 1
    
25.         elif x != a:
    
26.             return 0
    
27.     elif n + 1 >= 0:
    
28.         if x < a:
    
29.             return 0
    
30.         elif x >= a:
    
31.             return P_Moment * (x - a) ** ( n + 1 )        
    
32.         
    
33. # 引入用于奇异函数计算时所需的系数矩阵
    
34. 
    
35. # 创建外部载荷矩阵，数量、数值可根据具体案例修改
    
36. # 用于剪力计算的载荷系数
    
37. matrix_Load_For_ShearForce = np.array([73.75, -120, -45, -30, 30, 61.25])
    
38. # 用于弯矩计算的载荷系数，注意根据奇异函数规则表更新系数，与剪力的载荷系数有所不同
    
39. matrix_Load_For_Moment = np.array([73.75, -120, -45, -15, 15, 61.25])
    
40. 
    
41. # 创建外部载荷位置矩阵，数量、数值可根据具体案例修改
    
42. matrix_Load_Loaction = np.array([0, 2, 4, 6, 9, 12])
    
43. 
    
44. # 创建阶次矩阵，数量、数值可根据具体案例修改
    
45. matrix_n = np.array([0, -1, 0, 1, 1, 0])
    
46. 
    
47. # 定义横坐标x的起点、终点、细分数量,细分数量影响计算精度
    
48. x = np.linspace(0, 12, 5000)
    
49. # 为shear_force_values生成与x矩阵相同的全零数组
    
50. shear_force_values = np.zeros_like(x)
    
51. 
    
52. # 为moment_values生成与x矩阵相同的全零数组
    
53. moment_values = np.zeros_like(x)
    
54. 
    
55. # 开始遍历矩阵
    
56. for i in range(len(matrix_Load_Loaction)):
    
57.     P_Shear = matrix_Load_For_ShearForce[i] # 创建P并将负载矩阵中的序号为i元素依次赋值给P
    
58.     P_Moment = matrix_Load_For_Moment[i] # 创建P并将负载矩阵中的序号为i元素依次赋值给P
    
59.     a = matrix_Load_Loaction[i] # 创建a并将负载位置矩阵中的序号为i元素依次赋值给a
    
60.     n = matrix_n[i] # 创建n并将阶次矩阵中的序号为i元素依次赋值给n
    
61. 
    
62. # 将shear_force的函数变成一个可以接受数组作为输入的向量化函数shear_force_values  
    
63. # 将x、a和n作为参数传递给向量化的shear_force函数。
    
64. # 向量化函数将对每个元素进行计算，并返回一个具有相同形状的数组，其中包含计算后的结果
    
65. # shear_force_values +=将上述乘积添加到名为shear_force_values的变量中
    
66. # 如果shear_force_values之前已经有值，这将进行逐元素的相加。
    
67. # 这样，shear_force_values逐渐累积了每个循环迭代中的剪力值
    
68.     shear_force_values += np.vectorize(shear_force)(x, a, n)
    
69.    
    
70. # 同上
    
71.     moment_values += np.vectorize(moment_value)(x, a, n)
    
72. 
    
73. # 这行代码的目的是计算剪力值，并将结果累积在变量shear_force_values中
    
74. # 通过使用向量化函数和广播机制，可以高效地对输入数组的每个元素进行计算
    
75. # 并利用向量化的特性，避免了显式地编写循环来处理数组的每个元素
    
76.    
    
77. def plot_shear_force_diagram(x, shear_force_values):
    
78.     plt.plot(x, shear_force_values)
    
79.     plt.xlabel('Position (m)')
    
80.     plt.ylabel('Shear Sorce (kN)')
    
81.     plt.title('Shear Force Diagram')
    
82.     plt.grid(True)
    
83. 
    
84. def plot_moment_diagram(x, moment_values):
    
85.     plt.plot(x, moment_values)
    
86.     plt.xlabel('Position (m)')
    
87.     plt.ylabel('Moment Values (kN·m)')
    
88.     plt.title('Moment Diagram')
    
89.     plt.grid(True)
    
90.    
    
91. # 创建一个包含两个子图的图形
    
92. fig, axs = plt.subplots(2)
    
93. 
    
94. # 在第一个子图中绘制剪力图
    
95. plt.sca(axs[0])
    
96. plot_shear_force_diagram(x, shear_force_values)
    
97. plt.ylim(-100, 100)
    
98. plt.yticks(np.arange(0, 120, 40))
    
99. 
    
100. # 在第二个子图中绘制弯矩图
     
101. plt.sca(axs[1])
     
102. plot_moment_diagram(x, moment_values)
     
103. plt.ylim(-100, 100)
     
104. plt.yticks(np.arange(0, 500, 100))
     
105. 
     
106. # 调整子图之间的间距
     
107. plt.subplots_adjust(hspace=1)
     
108. 
     
109. # 显示图形
     
110. plt.grid(True)
     
111. plt.show()

复制代码

零件工程师

再优化一下老的Mathcad文件，总得来讲还是Mathcad编起来简单的多，然而Python免费。


哈哈，从Excel，到Mathcad，最终到Python，这个帖子算是完结了，问题结束。