简介：这份资源围绕「建筑地震保护系统」的建模与分析展开，面向建筑工程、结构动力学方向的科研人员及高年级本科生、研究生，适用于地震带新建建筑的抗震方案设计与技术预研。内容以弹性隔振层、榫头与卡槽自动锁定机构、弹簧-阻尼系统为核心，先建立系统在小振幅状态下的动力学模型，再讨论不同参数取值下榫头锁入卡槽的条件与锁入后的运动规律，并额外设计一种以扭转方式抵抗侧向风载的机械装置，给出原理图、模型方程与舒适性指标对设计参数的要求。压缩包仅含1个PDF文件，约289KB，为任务书与报告规范文档，附有系统结构示意图、评分标准及A4排版、Mathtype公式编辑等格式与纪律要求。已有69人学习下载，可帮助读者理解隔震与限位锁定的耦合建模思路，掌握从参数假设、解析求解到Matlab数值仿真的完整流程，并借鉴其图表规范与结论表达方式。
1. 从一次横波输入说起：这套隔震-锁榫系统到底在算什么
地基横波过来的时候，建筑很少是被“晃倒”的，多数破坏来自两件事：隔震层位移顶到限位后刚度突变，以及上部结构的楼层加速度峰值被拉爆。这套建筑地震保护系统给出的思路是分阶段处理——小振幅阶段让弹性隔震层按设计意图工作，把地面加速度隔开；振幅加大后，建筑与槽座的相对位移超过榫头与卡槽之间的间隙，榫头压入卡槽锁死，建筑和槽座连成一体，再由弹簧-阻尼这一对组合去耗能，用多出来的一个耦合自由度把振幅压下来。
真正要算清楚的是三个问题：小振幅时榫头会不会误锁（误锁等于隔震层失效，加速度直接倒灌进上部结构）、大振幅时锁不锁得上、锁上以后是更快衰减还是把位移顶得更大。适合动手做的人包括结构动力学课程设计的高年级学生、研究生，以及要给隔震层配限位与耗能装置的工程师。集中质量模型用 MATLAB 的 ode45 跑分段系统就够用，三维壳单元模型留着校核前几阶频率，不拿来做参数扫描。
2. 小振幅工况：二自由度隔震模型与 ode45 状态空间实现
2.1 从示意图拆自由度与参数符号
示意图里九个部件，进入运动方程的其实只有四个：地基输入、建筑集中质量、槽座支撑杆集中质量，以及两套弹簧-阻尼对。滚轮的作用是把槽座的运动约束成一维平动，接触摩擦先忽略；榫头和卡槽在小振幅阶段不接触，只贡献一个几何间隙 δ。自由度划错，后面分段积分的切换时刻就全错。
| 符号 | 含义 | 典型取值 | 单位 |
|---|---|---|---|
| m_b | 建筑集中质量 | 4.0e5 | kg |
| m_c | 槽座支撑杆质量 | 6.0e4 | kg |
| k_1、c_1 | 弹性隔震层刚度、阻尼 | T_b = 2.0 s 反算 | N/m、N·s/m |
| k_2、c_2 | 槽座弹簧、阻尼器参数 | T_c = 1.6 s 反算 | N/m、N·s/m |
| δ | 榫头与卡槽单侧间隙 | 0.02 ~ 0.12 | m |
| u_g、ü_g | 地基位移、绝对加速度 | 0.1g ~ 0.6g | m、m/s² |
用集中质量模型而不是直接上三维有限元，理由很实际：横波输入下隔震结构的响应由前两三阶模态主导，集中质量模型能解析地给出刚度、阻尼的物理含义，一轮扫参几秒钟跑完；ANSYS 或 SAP2000 的壳/梁单元模型只在最后校核频率和局部受力，误差控制在 5% 以内即可，不参与调参。
2.2 相对坐标下的运动方程
取相对地基的位移 x_b = u_b − u_g、x_c = u_c − u_g，两个子系统在小振幅阶段完全解耦：
m_b ẍ_b + c_1 ẋ_b + k_1 x_b = −m_b ü_g m_c ẍ_c + c_2 ẋ_c + k_2 x_c = −m_c ü_g
用相对坐标而不是绝对坐标，好处是数值上不必携带一个幅值远大于结构响应量的地面位移，ode45 的绝对误差容限可以放心收紧到 1e-11 量级。橡胶隔震层的滞回耗能用等效粘滞阻尼比折算，ζ_1 取 0.05 ~ 0.15；槽座那套弹簧-阻尼器是机械式的，ζ_2 可以取到 0.10 ~ 0.20。两边阻尼比不用凑成一样，后面的调谐判据会说明为什么。
2.3 状态空间与 ode45 实现
把两个二阶方程拆成四个一阶方程，写成匿名函数交给 ode45：
function dx = smallAmp(t, x, p) % x = [x_b; v_b; x_c; v_c]，四个量都相对地基 ag = -p.Ag * sin(p.w * t); % 地基绝对加速度，正弦横波输入 dx = zeros(4,1); dx(1) = x(2); dx(2) = (-p.mb*ag - p.c1*x(2) - p.k1*x(1)) / p.mb; dx(3) = x(4); dx(4) = (-p.mc*ag - p.c2*x(4) - p.k2*x(3)) / p.mc; end驱动脚本里参数一次配齐，单位统一到 kg、m、s、N：
p.mb = 4.0e5; p.mc = 6.0e4; p.k1 = p.mb*(2*pi/2.0)^2; % 隔震层，对应 T_b = 2.0 s p.k2 = p.mc*(2*pi/1.6)^2; % 槽座弹簧，对应 T_c = 1.6 s p.c1 = 2*0.08*sqrt(p.mb*p.k1); % ζ_1 = 0.08 p.c2 = 2*0.10*sqrt(p.mc*p.k2); % ζ_2 = 0.10 p.Ag = 0.15*9.81; p.w = 2*pi/1.6; p.delta = 0.05; [t, x] = ode45(@(t,x) smallAmp(t,x,p), [0 40], zeros(4,1)); drel = x(:,1) - x(:,3); % 榫头与卡槽的相对位移零初始条件表示从静止起振，40 s 足够走到稳态。-m·ag 这一项是地面加速度产生的惯性激励，符号别写反：正弦加速度取负号，对应位移从零向正向起摆。相对位移 drel 的峰值直接拿去和 δ 比，超了就说明小振幅假设不成立，必须切到锁入模型。
2.4 用传递函数预估锁入阈值与有限元校核
扫参之前先用解析式框一个范围，省掉大量盲跑。地面加速度到相对位移的传递函数是：
wb = sqrt(p.k1/p.mb); zb = p.c1/(2*sqrt(p.k1*p.mb)); wc = sqrt(p.k2/p.mc); zc = p.c2/(2*sqrt(p.k2*p.mc)); w = p.w; Hb = -1 / (wb^2 - w^2 + 2i*zb*wb*w); Hc = -1 / (wc^2 - w^2 + 2i*zc*wc*w); dX = abs(Hc - Hb) * p.Ag; % 稳态相对位移幅值估计 if dX < p.delta, disp('小振幅假设成立'); end这里藏着一个反直觉结论：当 wb = wc 且 zb = zc 时，Hb 与 Hc 完全相等，相对位移恒为零，榫头在任何振幅下都锁不上。工程上这既是好事也是坏事——想让榫头在大震时才动作，就把两个子系统的频率差当旋钮用，而不是一味加大 δ。算完解析值再回到有限元模型里取前两阶频率做对照，两边差 5% 以内就可以放心用集中质量模型扫参。
3. 榫头锁入卡槽的力学判据与事件检测
3.1 锁入需要同时满足的三个条件
只判“相对位移大于间隙”是不够的，回弹瞬间同样会穿过这个阈值，此时榫头没有压入的动量，硬判锁入会让仿真出现物理上不存在的刚度突变。完整的判据是三条同时成立：
| 条件 | 表达式 | 物理含义 | 不满足的后果 |
|---|---|---|---|
| 几何穿透 | |x_b − x_c| ≥ δ | 榫头够到卡槽 | 根本没接触 |
| 速度同向 | (ẋ_b − ẋ_c)·(x_b − x_c) > 0 | 正在继续压入，而非回弹 | 假锁，位移被凭空截断 |
| 接触受压 | F > 0 | 卡槽能推不能拉 | 立即脱开，出现反复切换 |
第三条在后处理里反算就行，前两条必须在事件函数里判。
3.2 用 odeset 的 Events 精确定位锁入时刻
让 ode45 自己找到穿越点，锁入时刻的精度直接取决于事件函数的写法：
opts = odeset('RelTol',1e-9, 'AbsTol',1e-11, 'Events', @lockEvent); function [val, ist, dir] = lockEvent(~, x, p) val = x(1) - x(3) - p.delta; % 正向穿越：相对位移继续增大 ist = 1; % 触发后停止积分，交给第二段接手 dir = 1; % 反向穿越另写一个事件，不要用 abs() end用 abs() 包住事件函数是常见的坑：绝对值在零点不可导，direction 参数失去意义，事件求解器会在零点附近反复缩短步长，步长被压到 1e-12 量级还触发不了，跑一次要好几分钟。正确做法是正向、反向各写一个事件函数，分别设 dir = 1 和 dir = −1。如果你采用近乎刚性的接触刚度模拟压入过程，记得把 Mass 矩阵显式传给 odeset，否则 ode45 会把这个问题当非刚性问题处理。
3.3 用小振幅段做一次能量核对
换地震波或者改参数之后，第一步不是看图，而是核对能量，确认数值耗散没有偷走响应：
Ein = -trapz(t, p.mb*ag.*x(:,2) + p.mc*ag.*x(:,4)); % 地面输入能量 Edis = trapz(t, p.c1*x(:,2).^2 + p.c2*x(:,4).^2); % 阻尼耗散能量 res = (Ein - Edis) / max(abs(Ein));ag 需要按时间向量重算一遍再相乘。res 应该落在 1e-3 以内；如果超过 1%，说明容限太松或者正弦频率与某个子系统频率撞上了，先把 RelTol 收到 1e-10 再确认一次。锁入前这一段是全流程里最干净的部分，这里对不上，后面分段模型的结果都不必看。
4. 锁入后的耦合模型：刚度切换与残余位移
4.1 约束方程与合并自由度
榫头锁进卡槽之后，两个质量之间只剩刚性约束：x_b = y + δ、x_c = y，y 是槽座相对地基的位移。把 2.2 节的两个方程相加，接触力 F 作为内力自动消掉：
(m_b + m_c) ÿ + (c_1 + c_2) ẏ + (k_1 + k_2) y = −(m_b + m_c) ü_g − k_1 δ
两个直接结论：锁入后系统退化成单自由度，质量取和、刚度和阻尼都取和；右端多出一项常数 −k_1 δ，相当于给系统加了一个静力偏置，稳态位置不在零而在 −k_1δ/(k_1+k_2)。这一项决定了残余位移，设计卡槽行程时要把它算进去。
4.2 锁入瞬间的初值与初加速度
第二段积分的初值来自第一段事件的终态，别用地面坐标去减：
y0 = [xe(end,3); xe(end,4)]; % 槽座位移与速度 [t2, y2] = ode45(@(t,y) locked(t,y,p), [te(end) 60], y0); function dy = locked(t, y, p) ag = -p.Ag*sin(p.w*t); dy = zeros(2,1); dy(1) = y(2); dy(2) = (-(p.mb+p.mc)*ag - (p.c1+p.c2)*y(2) ... - (p.k1+p.k2)*y(1) - p.k1*p.delta) / (p.mb+p.mc); end初加速度不用额外指定，ode45 会自己由方程算出来。需要确认的是切换前后建筑位移是否连续——理论上 x_b = y + δ 会有一个 δ 的跳跃，这个跳跃是物理的，来自榫头压入的瞬间；如果你在结果里看到建筑位移连续而槽座跳了，说明 δ 的符号取反了。
4.3 接触力反算与解锁判据
接触面只能推不能拉，反算 F 是判断会不会脱开的唯一手段：
ag2 = -p.Ag*sin(p.w*t2); acc = (-(p.mb+p.mc).*ag2 - (p.c1+p.c2).*y2(:,2) ... - (p.k1+p.k2).*y2(:,1) - p.k1*p.delta) / (p.mb+p.mc); F = p.mb*(acc + ag2) + p.c1.*y2(:,2) + p.k1.*(y2(:,1) + p.delta);F 是卡槽作用在建筑上的接触力，F ≤ 0 表示接触面被拉开，榫头脱出，此时应当回到 2 节的解耦模型重新起算，而不是继续按合并自由度积分。反复脱开-再锁在真实结构里意味着撞击，加速度谱会出现高频尖峰，这也是限制接触刚度不能取太小的原因。
4.4 加锁前后的关键指标对比
| 指标 | 小振幅解耦（δ = 0.05 m） | 锁入后合并 | 说明 |
|---|---|---|---|
| 系统自由度 | 2 | 1 | 刚度切换的直接结果 |
| 等效周期 | 2.0 s / 1.6 s 两个 | 由 k_1+k_2、m_b+m_c 决定，通常短于 1.6 s | 周期变短，加速度抬升 |
| 位移峰值 | 由相对位移包络控制 | 由静力偏置与瞬态叠加 | 看是否超过卡槽行程 |
| 稳态位置 | 零 | −k_1δ/(k_1+k_2) | 残余位移来源 |
| 衰减速度 | 各自按 ζ 衰减 | 按合并后的等效阻尼比衰减 | 通常快于隔震层单独作用 |
加速度抬升是锁入的代价，换来的是一次性把位移封顶。扫参时这两条曲线必须放在同一张图上看，只盯位移或者只盯加速度都会给出错误的推荐值。
5. 参数敏感性：间隙、刚度比与阻尼比的扫参实现
5.1 扫什么、看什么
参数里有三个真正可调的旋钮：榫头间隙 δ、槽座与隔震层的频率比（通过 k_2 体现）、两边的阻尼比。对应四个评价指标：是否锁入、锁入时刻 t_lock、建筑绝对加速度峰值 a_peak、末态残余位移。把两段积分封成一个函数，每轮返回这一个结构体：
function M = simulateTwoPhase(p) M = struct('lock',0,'tlock',NaN,'dmax',0,'apeak',0,'res',0,'Fmin',Inf); opts = odeset('RelTol',1e-9,'AbsTol',1e-11,'Events',@lockEvent); [t1,x1,te,xe] = ode45(@(t,x) smallAmp(t,x,p), [0 60], zeros(4,1), opts); M.dmax = max(abs(x1(:,1)-x1(:,3))); if isempty(te), return; end % 全程没锁上，直接返回 M.lock = 1; M.tlock = te(end); y0 = [xe(end,3); xe(end,4)]; [t2,y2] = ode45(@(t,y) locked(t,y,p), [te(end) 60], y0); % ……此处接 4.3 节反算 acc 与 F，取峰值与最小值 end60 s 的积分时长是留给稳态的余量；Fmin 用来判断稳定锁入还是反复脱开，比看位移曲线直观得多。
5.2 参数取值与扫参脚本
| 参数 | 扫描范围 | 依据 |
|---|---|---|
| δ | 0.02 ~ 0.12 m | 覆盖卡槽常见机械行程 |
| T_c | 1.0 ~ 2.5 s | 与隔震层周期形成频率差 |
| k_2/k_1 | 0.5 ~ 2.0 | 刚度比决定锁入后的等效周期 |
| ζ_1 / ζ_2 | 0.05 ~ 0.15 / 0.10 ~ 0.20 | 橡胶层与机械阻尼器的典型区间 |
| A_g | 0.1g ~ 0.6g | 小震到大震的输入幅值 |
dList = 0.02:0.02:0.12; rList = 0.5:0.25:2.0; for i = 1:numel(dList) for j = 1:numel(rList) p.delta = dList(i); p.k2 = rList(j) * p.mc * (2*pi/1.6)^2; % 按 T_c = 1.6 s 等比缩放 p.c2 = 2*0.10*sqrt(p.mc*p.k2); % 刚度变了，阻尼必须重算 T{i,j} = simulateTwoPhase(p); end endk_2 缩放时 c_2 一定要跟着重算，否则等效阻尼比会随刚度漂移，扫出来的趋势图没法横向比较，这是整个扫参里最容易被忽略的一步。
5.3 结果怎么读、哪里会翻车
- δ 偏小：小震就把榫头锁上，隔震层等于被旁路，a_peak 明显抬升，而且锁入时刻提前到激励上升段。这一侧的曲线特征是 a_peak 随 δ 减小而单调增大。
- δ 偏大：大震下相对位移够不到间隙，全程按解耦模型走，位移峰值靠隔震层自己的阻尼压不下来。中间存在一个窗口，窗口宽度由两个子系统的频率差决定。
- 刚度比接近 1 且阻尼匹配：相对位移被压到最小，误锁概率最低，但锁入后的合并刚度也最低，残余位移偏大，需要和卡槽行程一起折中。
- 阻尼取得过大：相对速度被压小，判据里的“速度同向”条件更容易在临界点附近反复触发，事件求解器会频繁中断。
- 数值上还要防两个坑：接触刚度取得过小会让接触力出现负值，模型误报脱开；事件方向设反会导致同一个穿越点被重复触发，积分卡在原地。
6. 扭转抗风装置与舒适度指标的校核方法
第（4）问要的是一套把横向风振转成扭转运动的装置。常见做法是用滚珠丝杠或齿条-齿轮把建筑顶部的水平摆动换成飞轮旋转，飞轮转动惯量 J 经传动比 i（单位 rad/m）折算成等效质量 m_eq = i²J，再配一根扭簧构成扭转调谐质量阻尼器：
m_s ẍ_s + c_s ẋ_s + k_s x_s = F_wind + c_t(θ̇ − i ẋ_s) + k_t(θ − i x_s) J θ̈ + c_t(θ̇ − i ẋ_s) + k_t(θ − i x_s) = 0
调谐初值按经典 TMD 公式给，调谐比 f 取 0.9 ~ 1.0，扭转阻尼比 ζ_t 取 0.05 ~ 0.15，惯容比 μ = m_eq/m_s 由加速度降幅要求反算：
ws = sqrt(p.ks/p.ms); f = 0.95; % 调谐比 zt = 0.08; p.kt = f^2 * ws^2 * p.J; % 扭簧刚度 N·m/rad p.ct = 2*zt*f*ws*p.J; % 扭转阻尼 N·m·s/rad mu = 0.03; p.i = sqrt(mu * p.ms / p.J); % 传动比，使 m_eq = i^2*J = mu*m_s舒适度指标用两个：10 分钟时程内的峰值加速度 a_peak 和均方根加速度 a_rms。住宅、公寓类建筑 a_peak 控制在 0.15 m/s² 以内，办公、旅馆类控制在 0.25 m/s² 以内，校核时程取 10 年重现期风荷载。
| 建筑类型 | a_peak 限值 (m/s²) | 校核时长 | 对设计参数的约束 |
|---|---|---|---|
| 住宅、公寓 | 0.15 | 10 min | μ 需更大或 ζ_t 提到 0.10 以上 |
| 办公、旅馆 | 0.25 | 10 min | μ = 0.02 ~ 0.03 通常够用 |
| 传动机构 | — | — | 飞轮转速 ω_max = i·v_max 须低于轴承许用转速 |
验证顺序建议先把解析传递率扫一遍频，确认峰值降幅达到 40% ~ 60% 且共振峰没有向低频平移过多，再用 Davenport 或 Kaimal 谱生成风时程做时域复核。最常翻车的地方不在调谐：飞轮许用转速是硬约束，先把 ω_max = i·v_max 算出来卡住传动比 i 的上限，再回头选惯量 J，往往比先挑飞轮尺寸再凑 μ 少走两轮返工。
本文还有配套的精品资源，点击获取
