利用Python代码对一些力学上的定义和概念进行可视化

想当年,我因为种种原因从理科基地班被人清退了,跑到了普通的班级。有一门课,我大半学期都没赶上,可期末还得参加考试,于是就去听了几次课。那老师讲课,口音实在太重了,把"挠度"说成"老陆",我基本上是听不懂的。当时还有个同学说他"嘴里有棉裤"。现在想想,这么说人家其实不太好,可当时确实是真听不明白。

说起来,现在有些高校在申请教师资格的时候,好像已经不要求普通话二级甲等以上了?有博士学位是不是就免了?记不太清这个标准了。但不管怎么说,还是希望高校教师能尽量把普通话说清楚,至少得让学生听懂。当时那位老师讲的东西,我好多都理解得不透,好在对付着考试算是过去了。

不过大学里应对考试是一回事,真学没学到东西又是另一回事。有的老师不负责任,直接给原题。有的学生觉得这老师好啊,给原题能过,然后就让大家背答案,ABCD对应1234,还真有孩子就背了。这可不是好事,到最后就是自己糊弄自己,根本没学明白,把1234背下来往那儿一放,既糊弄自己又糊弄家长。可孩子们那会儿不懂事,还觉得老师好。这种现象,除了抨击一下,也不好说别的。

咱只能想想自己当年没学好。后来我就打算再了解一下,补补课。这一补,就有了今天这些内容。本来想着严格对应地质现象来讲,可那样太专门、太枯燥,还显得有点太端着。何况我已经离地学很远了,所以就通用地讲讲这些概念,用代码画图,把力学里这些定义和概念直观地摆出来。

顺便说一句,这篇也想做个示范。我学的专业跟计算机离得挺远,可这些年靠着学的那点编程,画图、算数、验证想法,省了不少事。所以说,就算不是计算机相关专业,学一点编程相关内容,对学自己本专业的东西,往往也是能帮上忙的。

那这些概念到底是啥呢?一说材料力学、固体力学,很多人就头大,觉得全是公式。其实这几个最基础的概念,全是身边的事儿:拉一根橡皮筋,拧一个瓶盖,用扳手拧螺丝,掰一根钢筋。咱们今天就把应力、应变、力矩、力偶这几个词一个个拆开,配图配上代码,一次讲明白。

应力:单位面积上分摊到的力

先想一个画面:一根橡皮筋,你用手往两边拉。你觉得"手在用力",可实际上,这股力不是集中在某一个点上,而是摊在整个截面上。把"力"除以"截面积",就是应力(stress),记作 \(\sigma\)

\[\sigma = \frac{F}{A}\]

单位是帕斯卡(Pa),也就是牛每平方米(N/m²)。一个大气压大概是十万帕,工程上常用兆帕(MPa)。

应力有两种最基础的:

  • 正应力 \(\sigma\):力垂直作用在面上,拉是拉应力,压是压应力。
  • 切应力 \(\tau\):力贴着面(相切)蹭过去,就像拿刀片刮表面那样。

这概念一点都不玄。你站着不动,体重压在两脚上,脚底承受的那个"压强",本质上就是一种正应力。穿细跟的鞋踩到别人脚,受力面积小,压强大——这就是应力大,挨踩的疼啊。

Image

这张图左边是正应力:一个正方形小块,左右两边受拉力往外拽;右边是切应力:上下两个面被"错着刮"。中间那个 \(\sigma\)\(\tau\),就是单位面积上分到的力。

下面的代码可以自己跑,画出同样的一张图:

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import os

plt.rcParams['font.sans-serif'] = ['Arial Unicode MS', 'Hiragino Sans GB', 'Heiti TC', 'Songti SC', 'Microsoft YaHei']
plt.rcParams['axes.unicode_minus'] = False
os.makedirs('images', exist_ok=True)

DARK = '#2c3e50'; RED = '#e74c3c'

def force_arrow(ax, x, y, dx, dy, color=RED, lw=2.2):
    ax.add_patch(mpatches.FancyArrowPatch(
        (x, y), (x+dx, y+dy), arrowstyle='-|>', mutation_scale=20, lw=lw, color=color))

fig, axes = plt.subplots(1, 2, figsize=(11, 5))
s = 1.0

ax = axes[0]  # 正应力
ax.add_patch(mpatches.Rectangle((-s, -s), 2*s, 2*s, fc='#dff0f7', ec=DARK, lw=2))
for y in (0.45, 0.0, -0.45):
    force_arrow(ax, -s, y, -1.0, 0)
    force_arrow(ax,  s, y,  1.0, 0)
ax.text(0, 0, 'σ', ha='center', va='center', fontsize=24, color=DARK)
ax.text(0, -2.0, '拉力 ÷ 截面积', ha='center', fontsize=13, color=RED)
ax.set_xlim(-3, 3); ax.set_ylim(-2.5, 2.5); ax.set_aspect('equal'); ax.axis('off')
ax.set_title('正应力 σ:力垂直作用在面上', fontsize=15)

ax = axes[1]  # 切应力
ax.add_patch(mpatches.Rectangle((-s, -s), 2*s, 2*s, fc='#fdf0d6', ec=DARK, lw=2))
for x in (-0.45, 0.0, 0.45):
    force_arrow(ax, x,  s, 0,  1.0)
    force_arrow(ax, x, -s, 0, -1.0)
ax.text(0, 0, 'τ', ha='center', va='center', fontsize=24, color=DARK)
ax.text(0, -2.0, '切向力 ÷ 截面积', ha='center', fontsize=13, color=RED)
ax.set_xlim(-3, 3); ax.set_ylim(-2.5, 2.5); ax.set_aspect('equal'); ax.axis('off')
ax.set_title('切应力 τ:力平行作用在面上', fontsize=15)

fig.suptitle('应力 = 单位面积上分摊到的力', fontsize=17, y=1.02)
fig.tight_layout()
plt.savefig('images/stress.png', dpi=150, bbox_inches='tight')

应变:变形量除以原来的尺寸

光有力还不够,还得看东西被"抻"了多少。一根橡皮筋原长 \(L\),拉长了 \(\Delta L\),那"变形程度"就是应变(strain):

\[\varepsilon = \frac{\Delta L}{L}\]

关键是除以原来的长度,这样不管东西本身多长,都能拿同一把尺子比较变形的剧烈程度。应变没有单位,是个比值。

跟应力对应,应变也有正应变和切应变:

  • 正应变 \(\varepsilon\):拉长或压短的比例。
  • 切应变 \(\gamma\):形状被"剪切"倾斜的程度,就是那个角度的变化量。

Image

左图:虚线是原来的样子,实线是拉长之后的,标了原长 \(L\) 和伸长量 \(\Delta L\)。右图:虚线是原来的方块,实线是被"错开"了的平行四边形,顶边整体滑移了一段,那个绿色的小角度 \(\gamma\) 就是切应变。

代码也不长:

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import os

plt.rcParams['font.sans-serif'] = ['Arial Unicode MS', 'Hiragino Sans GB', 'Heiti TC', 'Songti SC', 'Microsoft YaHei']
plt.rcParams['axes.unicode_minus'] = False
os.makedirs('images', exist_ok=True)

DARK = '#2c3e50'; RED = '#e74c3c'; BLUE = '#2980b9'; GRAY = '#95a5a6'; GREEN = '#2ecc71'

def force_arrow(ax, x, y, dx, dy, color=RED, lw=2.2):
    ax.add_patch(mpatches.FancyArrowPatch(
        (x, y), (x+dx, y+dy), arrowstyle='-|>', mutation_scale=20, lw=lw, color=color))

fig, axes = plt.subplots(1, 2, figsize=(11, 5))
L = 2.0; dL = 0.8

ax = axes[0]  # 正应变
ax.add_patch(mpatches.Rectangle((-L/2, -0.6), L, 1.2, fill=False, ls='--', ec=GRAY, lw=1.8))
ax.add_patch(mpatches.Rectangle((-L/2, -0.6), L+dL, 1.2, fc='#dff0f7', ec=BLUE, lw=2))
for y in (0.3, -0.3):
    force_arrow(ax, -L/2, y, -0.9, 0)
    force_arrow(ax,  L/2+dL, y, 0.9, 0)
ax.annotate('', xy=(L/2+dL, 1.35), xytext=(L/2, 1.35),
            arrowprops=dict(arrowstyle='<->', color=BLUE, lw=1.6))
ax.text(L/2+dL/2, 1.6, 'ΔL(伸长量)', ha='center', fontsize=12, color=BLUE)
ax.annotate('', xy=(0.9, -1.1), xytext=(-0.9, -1.1),
            arrowprops=dict(arrowstyle='<->', color=GRAY, lw=1.6))
ax.text(0, -1.5, '原长 L', ha='center', fontsize=12, color=GRAY)
ax.text(L/2+dL+0.3, 0, 'ε = ΔL / L', ha='left', va='center', fontsize=14, color=DARK)
ax.set_xlim(-3.3, 5.5); ax.set_ylim(-2.2, 2.2); ax.set_aspect('equal'); ax.axis('off')
ax.set_title('正应变 ε:拉长(或压短)的比例', fontsize=15)

ax = axes[1]  # 切应变
shift = 0.9
ax.add_patch(mpatches.Rectangle((-L/2, -0.6), L, 1.2, fill=False, ls='--', ec=GRAY, lw=1.8))
corners = [(-L/2, -0.6), (L/2, -0.6), (L/2+shift, 0.6), (-L/2+shift, 0.6)]
ax.add_patch(mpatches.Polygon(corners, closed=True, fc='#fdf0d6', ec=BLUE, lw=2))
for x in (0.0, 0.45, 0.9):
    force_arrow(ax, x-0.6, 0.6, 0, 0.85)
    force_arrow(ax, x-0.6+L, -0.6, 0, -0.85)
ang = np.degrees(np.arctan2(shift, L))
ax.add_patch(mpatches.Arc((-L/2, -0.6), 0.9, 0.9, theta1=0, theta2=ang, color=GREEN, lw=2))
ax.text(-L/2+0.55, -0.6+0.22, 'γ', fontsize=16, color=GREEN)
ax.text(-L/2+shift/2, 1.35, '上边沿横向滑移', ha='center', fontsize=12, color=BLUE)
ax.text(L/2+1.9, 0, 'γ = 滑移量 / 高度', ha='left', va='center', fontsize=13, color=DARK)
ax.set_xlim(-2.2, 5.8); ax.set_ylim(-2.2, 2.4); ax.set_aspect('equal'); ax.axis('off')
ax.set_title('切应变 γ:形状发生"剪切"倾斜', fontsize=15)

fig.suptitle('应变 = 变形量与原来尺寸的比值', fontsize=17, y=1.02)
fig.tight_layout()
plt.savefig('images/strain.png', dpi=150, bbox_inches='tight')

应力除以应变,就是弹性模量 \(E\)\(E = \sigma / \varepsilon\)。它衡量材料"有多硬"——同样抻它,应力涨得快的,就是难拉动的硬材料。在弹性阶段,胡克定律成立:\(\sigma = E\varepsilon\)。这不就是中学学的弹簧嘛,F 换成 σ,伸长量换成 ε,比例系数换成 E。

应力-应变曲线:每种材料各有各的"脾气"

把一根材料样品慢慢拉,一边记录应力一边记录应变,画出来的就是应力-应变曲线。这条曲线基本就是材料的"性格画像"。

Image

图里粗略画了三种典型材料:

  • 低碳钢(韧性材料):先是笔直的弹性段,撤力就复原;然后到屈服点 \(\sigma_y\),过了这点材料就"回不去了",产生永久变形,曲线上有个平台;再往后是强化阶段,能扛的应力继续涨到抗拉强度 \(\sigma_u\);最后某处开始变细(颈缩),没扛多久就断了。整个过程"能屈能伸",断之前有大量变形预警,这就是韧性。
  • 铸铁(脆性材料):几乎是一条直线,没怎么变形,"咔嚓"一下就断了。看着硬,其实很脆,没有预警。
  • 橡胶(高弹性材料):弹性模量很小,软,能抻出好几倍的长度不断,曲线是条往上翘的弧线。
import numpy as np
import matplotlib.pyplot as plt
import os

plt.rcParams['font.sans-serif'] = ['Arial Unicode MS', 'Hiragino Sans GB', 'Heiti TC', 'Songti SC', 'Microsoft YaHei']
plt.rcParams['axes.unicode_minus'] = False
os.makedirs('images', exist_ok=True)

DARK = '#2c3e50'; BLUE = '#2980b9'; RED = '#e74c3c'; GREEN = '#2ecc71'
fig, ax = plt.subplots(figsize=(11, 6))

# 低碳钢:弹性→屈服平台→强化→颈缩→断裂
ey, sy = 0.002, 400.0
E = sy / ey
e1 = np.linspace(0, ey, 60); s1 = E * e1
e2 = np.linspace(ey, 0.02, 40); s2 = np.full_like(e2, sy)
t = np.linspace(0, 1, 120) ** 0.9
e3 = 0.02 + (0.18 - 0.02) * t; s3 = sy + (600 - sy) * t
t2 = np.linspace(0, 1, 60) ** 1.4
e4 = 0.18 + (0.22 - 0.18) * t2; s4 = 600 - (600 - 500) * t2
ax.plot(np.concatenate([e1, e2, e3, e4]), np.concatenate([s1, s2, s3, s4]),
        color=BLUE, lw=2.4, label='低碳钢(韧性材料)')

# 铸铁(脆性)
e_ci = np.linspace(0, 0.0045, 60)
s_ci = np.clip(38000 * e_ci * (1 - 8 * e_ci), 0, 140)
ax.plot(e_ci, s_ci, color=RED, lw=2.4, label='铸铁(脆性材料)')

# 橡胶(高弹性)
e_r = np.linspace(0, 4.0, 120)
ax.plot(e_r, 3.2 * (np.exp(0.75 * e_r) - 1), color=GREEN, lw=2.4, label='橡胶(高弹性材料)')

ax.scatter([ey], [sy], color=BLUE, s=40, zorder=5)
ax.annotate('屈服点 σy', xy=(ey, sy), xytext=(0.05, 620),
            arrowprops=dict(arrowstyle='->', color=DARK), fontsize=12)
ax.annotate('抗拉强度 σu', xy=(0.18, 600), xytext=(0.13, 500),
            arrowprops=dict(arrowstyle='->', color=DARK), fontsize=12)
ax.annotate('断裂', xy=(0.22, 500), xytext=(0.13, 180),
            arrowprops=dict(arrowstyle='->', color=DARK), fontsize=12)

ax.text(0.0004, 300, '弹性阶段\n(撤力就复原)', fontsize=11, color=BLUE)
ax.text(0.10, 520, '强化阶段', fontsize=12, color=BLUE)
ax.text(0.20, 360, '颈缩', fontsize=12, color=BLUE)
ax.text(0.0004, 100, '脆性材料:\n没怎么变形就断了', fontsize=11, color=RED)
ax.text(1.6, 40, '橡胶:\n能抻好几倍不断', fontsize=11, color=GREEN)

ax.set_xlabel('应变 ε', fontsize=13)
ax.set_ylabel('应力 σ(MPa)', fontsize=13)
ax.set_title('应力-应变曲线:不同材料各有各的"脾气"', fontsize=15)
ax.set_xlim(0, 4.3); ax.set_ylim(0, 700); ax.grid(True, alpha=0.3)
ax.legend(loc='lower right', fontsize=11)
fig.tight_layout()
plt.savefig('images/ss_curve.png', dpi=150, bbox_inches='tight')

工程上选材料,就是在看这条曲线:要扛得住别断,看强度;要能变形缓冲不断,看韧性;要硬,看弹性模量。

力矩:让物体转动的"扭转本领"

前面说的力,是让物体平移的。可拧螺丝、开门、转方向盘,需要的不是"拽",而是"拧"。拧的本事叫力矩(moment,也叫扭矩),公式是:

\[M = F \times d\]

\(F\) 是力,\(d\)力臂——从转动轴到力的作用线的垂直距离。单位是牛米(N·m)。

Image

图上这个扳手最好理解了:螺母转轴在左端,你在扳手末端用力 F 往下压。力臂 d 就是转轴到那条竖直力作用线的垂直距离。力矩让扳手绕着螺母转。

这里面有个特别实用的结论:力臂越长越省力。所以加长扳手(或者套一根管子加长力臂)就能撬动更紧的螺丝;门把手都装在离合页远远的一边,开起来才省劲。这些都是力矩 \(M = F \times d\) 在起作用。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import os

plt.rcParams['font.sans-serif'] = ['Arial Unicode MS', 'Hiragino Sans GB', 'Heiti TC', 'Songti SC', 'Microsoft YaHei']
plt.rcParams['axes.unicode_minus'] = False
os.makedirs('images', exist_ok=True)

DARK = '#2c3e50'; RED = '#e74c3c'; BLUE = '#2980b9'; GREEN = '#2ecc71'; GRAY = '#95a5a6'
fig, ax = plt.subplots(figsize=(10, 6))

ax.add_patch(mpatches.Circle((0, 0), 0.32, fc='#bdc3c7', ec=DARK, lw=2))
ax.add_patch(mpatches.Circle((0, 0), 0.14, fc=DARK))
ax.plot([0, 3.4], [0, 0.65], color=DARK, lw=8, solid_capstyle='round')

fx, fy = 3.4, 0.65
ax.add_patch(mpatches.FancyArrowPatch((fx, fy+0.9), (fx, fy-0.6),
             arrowstyle='-|>', mutation_scale=20, lw=3, color=RED))
ax.text(fx+0.18, fy+0.45, '力 F', fontsize=15, color=RED)
ax.plot([0, 3.4], [0, 0], color=BLUE, lw=1.4, ls='--')
ax.annotate('', xy=(3.4, 0.42), xytext=(3.4, -0.42),
            arrowprops=dict(arrowstyle='<->', color=BLUE, lw=1.6))
ax.text(3.5, 0.05, '力臂 d', fontsize=14, color=BLUE, ha='left')
# 转动方向:力向下压,扳手绕螺母顺时针转
th = np.linspace(np.radians(15), np.radians(-75), 80)
ax.plot(1.0*np.cos(th), 1.0*np.sin(th), color=GREEN, lw=2.2, solid_capstyle='round')
t = np.radians(-75)
ax.add_patch(mpatches.RegularPolygon((1.0*np.cos(t), 1.0*np.sin(t)), 3, radius=0.10,
             orientation=np.radians(np.degrees(np.arctan2(-np.cos(t), np.sin(t)))), color=GREEN))
ax.text(0.25, -0.75, 'M = F × d', fontsize=16, color=GREEN, fontweight='bold')
ax.text(0.3, -1.15, '(力臂越长越省力)', fontsize=12, color=GREEN)
ax.text(0, -2.0, '固定螺母', ha='center', fontsize=12, color=GRAY)
ax.set_xlim(-1.3, 5.2); ax.set_ylim(-2.4, 2.0); ax.set_aspect('equal'); ax.axis('off')
ax.set_title('力矩(扭矩):让物体转动的"扭转本领"', fontsize=15)
fig.tight_layout()
plt.savefig('images/moment.png', dpi=150, bbox_inches='tight')

力偶:一对力,只转不挪

有时候你会碰上一种特殊的"拧法":两个大小相等、方向相反的力,平行但不共线,一起作用在一个物体上。这一对力叫力偶(couple)。

Image

图上这块东西,左边一个力往上拽,右边一个大小一样、方向相反的力往下压。两股力加起来合力正好是零——所以它不会整体移动;但这一拽一压,恰好让它一个劲地转。这种"只转不挪"的效果,叫力偶矩:

\[M = F \times d\]

\(d\) 是那两个力之间的垂直距离。有意思的是,力偶矩跟你的参考点选在哪无关,只要这一对力在,拧的效果就是固定的。拧瓶盖、双手搓方向盘、拧毛巾,靠的全是力偶。

下面这张动图更直观,力偶作用在一个方块上,方块就这么转起来了,可它的位置一直没挪:

Image

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from matplotlib.animation import FuncAnimation
import os

plt.rcParams['font.sans-serif'] = ['Arial Unicode MS', 'Hiragino Sans GB', 'Heiti TC', 'Songti SC', 'Microsoft YaHei']
plt.rcParams['axes.unicode_minus'] = False
os.makedirs('images', exist_ok=True)

DARK = '#2c3e50'; RED = '#e74c3c'
fig, ax = plt.subplots(figsize=(5.0, 5.0))
frames = 30; W, H = 2.2, 1.0

def draw_frame(i):
    ax.clear()
    ang = -np.radians(i * 360 / frames)   # 左手上推、右手下压,整体顺时针转
    R = np.array([[np.cos(ang), -np.sin(ang)], [np.sin(ang), np.cos(ang)]])
    corners = np.array([[-W/2, -H/2], [W/2, -H/2], [W/2, H/2], [-W/2, H/2]])
    ax.add_patch(mpatches.Polygon(corners @ R.T, closed=True, fc='#dff0f7', ec=DARK, lw=2))
    p1 = np.array([-W/2, 0]) @ R.T
    p2 = np.array([ W/2, 0]) @ R.T
    v = np.array([0, 1]) @ R.T
    ax.add_patch(mpatches.FancyArrowPatch(p1, p1 + v, arrowstyle='-|>', mutation_scale=22, lw=2.6, color=RED))
    ax.add_patch(mpatches.FancyArrowPatch(p2, p2 - v, arrowstyle='-|>', mutation_scale=22, lw=2.6, color=RED))
    ax.text(0, 0, '力偶矩 M\n只转动,不移动', ha='center', va='center', fontsize=12, color=DARK)
    ax.set_xlim(-2.3, 2.3); ax.set_ylim(-2.3, 2.3); ax.set_aspect('equal'); ax.axis('off')

anim = FuncAnimation(fig, draw_frame, frames=frames, interval=60)
anim.save('images/couple_rotate.gif', writer='pillow', fps=15)

串起来记

把这几个概念用一句话串起来:

应力是单位面积上分摊到的力,应变是变形量相对于原尺寸的比例;力矩是力乘力臂的"拧"的本领,力偶是一对力只转不挪的拧法。

公式不用死记,把图看明白就行:应力应变说的是一块材料内部"扛了多少、变了多少",力矩力偶说的是"拿什么拧、怎么拧"。这几个概念弄明白了,回头再看受力、变形的各种例子,就不至于两眼一抹黑。

写这篇,除了讲清楚这几个概念,其实更想做个示范:就算不是计算机相关专业,也完全可以考虑学一点编程相关内容。拿它当个趁手的工具,画个图、算个数、验证一下想法,学自己专业的东西,往往也能帮上不少忙。这活儿不挑专业,门槛也没想象中那么高,有兴趣就能上手。

配图生成的完整代码就在上面的各个代码块里,复制就能跑。要是有讲得不清楚的地方,欢迎指正。