CN106909755A - 一种车辆荷载下沉管隧道管节竖向位移的计算方法 - Google Patents

Abstract

本发明提供一种车辆荷载下沉管隧道管节竖向位移的计算方法，采用Timoshenko梁来模拟管节，将地基等效为一系列并联的弹簧元件和阻尼元件，建立管节‑接头模型，接头模型中的抗剪单元和抗弯单元均由弹簧和阻尼并联组成；简化管节边界条件，将其考虑为自由‑自由；接头作用通过在管节端部添加集中力和集中弯矩实现，相邻端面所受集中力和集中弯矩大小相等、方向相反；本发明通过管节振型函数求解以及管节动力方程建立及求解，采用Newmark逐步积分法求解，得到第j段管节第n阶时间系数，结合管节模态函数能够得到管节纵向任意位置的竖向位移响应。

Description

一种车辆荷载下沉管隧道管节竖向位移的计算方法

技术领域

本发明属于地下工程技术领域，具体涉及一种车辆荷载下沉管隧道管节竖向位移的计算方法。

背景技术

沉管隧道建造于水底，需严格保证密封，因此对其管节端部位移控制有较高要求。然而，根据国内外工程案例，部分沉管隧道管节在运营期发生开裂并渗漏，对隧道正常使用带来极大危害。目前，国内外已建和在建沉管隧道多为公路隧道，因此，车辆荷载对管节位移影响较大。为探究车辆荷载对管节及接头位移响应的影响规律，保证隧道内行车的舒适性，同时预防管节接头发生渗漏，本发明将建立车辆荷载下沉管隧道竖向位移计算模型及方法。

目前，在国内外沉管隧道结构设计中大多将车辆荷载进行拟静力计算，有关车辆荷载对沉管隧道位移响应研究较少，且已有沉管隧道结构受力计算模型中，多将管节视为Euler梁，而未考虑其剪切变形。

发明内容

本发明的目的是克服现有技术中不足，提供一种车辆荷载下沉管隧道管节竖向位移的计算方法。

为了达到上述目的，本发明是通过以下技术方案实现的：

本发明考虑采用抗剪单元与抗弯单元并联来模拟柔性接头，具体地：本发明提供一种车辆荷载下沉管隧道管节竖向位移的计算方法，采用Timoshenko梁来模拟管节，将地基等效为一系列并联的弹簧元件和阻尼元件，建立管节-接头模型，如图1所示；接头模型中的抗剪单元和抗弯单元均由弹簧和阻尼并联组成；

简化管节边界条件，将其考虑为自由-自由；接头作用通过在管节端部添加集中力和集中弯矩实现，相邻端面所受集中力和集中弯矩大小相等、方向相反(如图2所示)；

本发明具体包括如下步骤：

步骤1)：管节振型函数求解

建立管节自由振动控制方程：

式中：κ为管节剪切系数，无量纲；

A为管节截面面积，单位为m²；

G为管节剪切模量，单位为Pa；

v为管节竖向位移，单位为m；

φ为管节转角，单位为rad；

ρ为管节密度，单位为kg/m³；

E为管节弹性模量，单位为Pa；

I为管节惯性矩，单位为m⁴；

x为距离管节端部的长度，单位为m；

t为时间，单位为s；

采用模态叠加法，假定管节竖向位移及转角表达式为：

式中：n为管节振动模态，无量纲；

ω_n为管节弯曲振动固有频率，单位为rad/s；

i为虚数单位；

me为所取最高管节模态数，无量纲；

将(2)代入(1)，并进行正交化解耦，令整理得到：

求解上述方程得到ω_n和λ_n(λ_1n,λ_2n)之间的关系：

将ω_n和λ_n(λ_1n,λ_2n)之间的关系代入位移v_n(x)和转角φ_n(x)的标准模态函数得到：

v_n(x)＝c_1nch(λ_1nx)+s_1nsh(λ_1nx)+c_2ncos(λ_2nx)+s_2nsin(λ_2nx) (5)

φ_n(x)＝c_1ng_1nsh(λ_1nx)+s_1ng_1nch(λ_1nx)-c_2ng_2nsin(λ_2nx)+s_2ng_2ncos(λ_2nx) (6)

式中：

c_1n、c_2n、s_1n、s_2n为振型函数系数；

根据管节简化模型建立边界条件：

式中：l为管节长度；

满足模态函数系数c_1n、c_2n、s_1n、s_2n不同时等于0，求解管节振动固有频率ω_n，从而得到管节模态振型，具体采用Matlab编程求解；上述方法适用于弹性体模态求解，而根据相关研究^[17]，自由边界条件下Timoshenko梁前两阶模态为刚体模态，其模态函数及频率为：

步骤2)：管节动力方程建立及求解

先建立管节受迫振动控制方程：

式中：F(x,t)为管节所受外力，单位N/m；

M(x,t)为管节所受弯矩N·m/m；

采用模态叠加法，假定梁的竖向位移及转角表达式为：

式中：q_n(t)为时间系数，单位为s；

将(10)代入(9)，进行正交化解耦得到第j段管节的第n阶振动常微分方程为：

式中：l_j为第j段管节长度，单位为m；

由于车辆质量相对管节质量可忽略不计，本文将车辆前后轴荷载等效为两个点源移动恒载：

P(t)＝∑P_mδ(x-(ut+x_m))δ(y) (12)

式中：P_m为时刻第m辆车作用在管节上的点荷载，单位为N；

δ(·)为狄拉克函数；

u为车辆行驶速度，单位为m/s；

x_m为车辆初始位置，单位为m；

y为管节横向坐标，单位为m；

假设车辆沿隧道轴线方向行驶，考虑车辆荷载、地基反力和接头集中力及弯矩作用，Fⁱ(x,t)和Mⁱ(x,t)的具体表达式为：

式中：k_j为接头抗剪单元弹簧系数，单位为N/m；

c_j为接头抗剪单元阻尼系数N·s/m；

k为地基等效弹簧系数，N/m²；

c为地基阻尼系数N·s/m²；

P_my为车辆等效横向均布荷载，单位为N/m；

式中：k_w为接头抗弯单元弹簧系数，单位为N·m/rad；

c_w为接头抗弯单元阻尼系数，单位为N·m·s/rad；

将(13)和(14)代入(11)最终得到：

将(15)整理成矩阵方程组，采用Newmark逐步积分法求解，得到第j段管节第n阶时间系数结合管节模态函数能够得到管节纵向任意位置的竖向位移响应。

与现有技术相比，本发明的有益效果如下：

现有参考文献的沉管隧道中大多将车辆荷载进行拟静力计算，有关车辆荷载引起管节竖向位移响应研究只有文献[1]有所涉及，但该研究只对单一管节在两端简支的情况下动力响应进行分析，而未考虑接头的影响。实际上，运营期间隧道管节端部存在竖向位移，且相邻管节间存在相互作用，因此，文献[1]提出的计算模型并不合理。本发明的理论基础扎实，考虑柔性接头对管节振动响应的影响，建立接头模型，同时考虑管节的弯曲变形和剪切变形，采用Timoshenko梁模拟管节，分析车辆荷载下管节竖向位移响应情况。实际计算中可利用Matlab软件编写程序，赋予车辆、管节结构、接头和地基合理的特性参数，借助计算机强大的运算能力提高计算速度和精度。

利用本发明提出的管节-接头竖向位移计算模型及计算方法，可对纵向坡度较小的管节车辆振动响应进行计算，从而研究管节及接头位移响应规律。此外，改变诸如车距、车速、车重、地基系数及接头系数等参数取值，计算不同工况下的管节动力响应结果并进行对比分析，可以研究单因素对管节的影响。

附图说明

图1为本发明中管节-接头模型示意图；

图2为本发明中管节纵向计算简化模型示意图；

图3为本发明中管节南端竖向位移示意图；

图4为本发明中管节北端竖向位移示意图；

图5为本发明中管节中点竖向位移示意图；

图6为本发明中接头两端竖向位移差示意图。

具体实施方式

下面结合说明书附图对本发明的技术方案作进一步说明：

如图1～6所示，本发明提出一种车辆荷载下沉管隧道管节竖向位移的计算方法的具体实施例，本发明以宁波甬江沉管隧道工程为背景，管节结构参数详见文献[2]。计算中取管节重度γ＝25kN/m³，弹性模量为E＝3.45×10⁴MPa，剪切系数κ＝π²/12^[19]，剪切模量G＝1.25×10⁴Mpa；取地基弹性抗力系数为k＝1×10⁴kN/m²，地基阻尼系数为c＝5×10²kPa·s。取接头抗剪刚度k_j＝1×10⁶kN/m，抗弯刚度k_w＝3×10⁶kN·m/rad，接头抗剪阻尼系数c_j＝5×10²kN·s/m，抗弯阻尼系数c_w＝5×10²kN·m·s/rad。假设车辆前后轴重P＝250kN，移动速度取隧道设计车速u＝60km/h，轴距l_w＝4m。车辆自北向南行驶，在隧道内均匀分布，相邻车辆间距取l＝13m。利用Matlab编写程序，求解车辆荷载下管节位移响应。由于甬江隧道管节模型分布基本对称，故取E3～E5管节和J4～J6接头两端位移差进行计算，结果如图3～6所示。管节位移以向下为正，向上为负。接头位移差以北端位移大于南端位移为正，反之为负。

图3～5结果表明，车辆荷载引起E3～5管节南北两端向上最大位移为1.4～2.6mm，引起E3～5管节中点向下最大位移为2.6～4.2mm。因此，车辆荷载引起管节中点位移幅度较端部位移幅度大，且越靠近江中段的管节振幅越大。图6结果表明，越靠近岸边的管节接头位移差受车辆荷载影响越大，最大位移差出现在J6接头，大小约为1.5mm。

上述实施例是对本发明的说明，不是对本发明的限定，任何对本发明简单变换后的方案均属于本发明的保护范围。

[1]苏勤卫.海底沉管隧道管段沉降与应变研究[D].杭州:浙江大学,2015.

[2]谢雄耀,王培,李永盛,等.甬江沉管隧道长期沉降监测数据及有限元分析[J].岩土力学,2014,35(8):2314-2324.

Claims (1)

1.一种车辆荷载下沉管隧道管节竖向位移的计算方法，其特征在于，采用Timoshenko梁来模拟管节，将地基等效为一系列并联的弹簧元件和阻尼元件，建立管节-接头模型，接头模型中的抗剪单元和抗弯单元均由弹簧和阻尼并联组成；

简化管节边界条件，将其考虑为自由-自由；接头作用通过在管节端部添加集中力和集中弯矩实现，相邻端面所受集中力和集中弯矩大小相等、方向相反；具体包括如下步骤：

步骤1)：管节振型函数求解

建立管节自由振动控制方程：

$\left\{\begin{array}{c}\mathrm{\κ}AG(\frac{{\∂}^{2}v(x,t)}{\∂x^2}-\frac{\∂\mathrm{\φ}(x,t)}{\∂x})=\mathrm{\ρ}A\frac{{\∂}^{2}v(x,t)}{\∂t^2}\\ EI\frac{{\∂}^2\mathrm{\φ}(x,t)}{\∂x^2}+\mathrm{\κ}AG(\frac{\∂v(x,t)}{\∂x}-\mathrm{\φ}(x,t\left)\right)=\mathrm{\ρ}I\frac{{\∂}^2\mathrm{\φ}(x,t)}{\∂t^2}\end{array}\right.---\left(1\right)$

式中：κ为管节剪切系数，无量纲；

A为管节截面面积，单位为m²；

G为管节剪切模量，单位为Pa；

v为管节竖向位移，单位为m；

φ为管节转角，单位为rad；

ρ为管节密度，单位为kg/m³；

E为管节弹性模量，单位为Pa；

I为管节惯性矩，单位为m⁴；

x为距离管节端部的长度，单位为m；

t为时间，单位为s；

采用模态叠加法，假定管节竖向位移及转角表达式为：

$\left\{\begin{array}{c}v\left(x,t\right)=\underset{n=1}{\overset{me}{\mathrm{\Σ}}}{v}_n\left(x\right)e^{{\mathrm{i\ω}}_{n}t}\\ \mathrm{\φ}\left(x,t\right)=\underset{n=1}{\overset{me}{\mathrm{\Σ}}}{\mathrm{\φ}}_n\left(x\right)e^{{\mathrm{i\ω}}_{n}t}\end{array}\right.---\left(2\right)$

式中：n为管节振动模态，无量纲；

ω_n为管节弯曲振动固有频率，单位为rad/s；

i为虚数单位；

me为所取最高管节模态数，无量纲；

将(2)代入(1)，并进行正交化解耦，令整理得到：

${{\mathrm{EI\λ}}_n}^4+{{\mathrm{\ω}}_n}^4(\frac{EI\mathrm{\ρ}}{\mathrm{\κ}G}+\mathrm{\ρ}I){{\mathrm{\λ}}_n}^2+({{\mathrm{\ω}}_n}^4\frac{{\mathrm{\ρ}}^{2}I}{\mathrm{\κ}G}-{{\mathrm{\ω}}_n}^2\mathrm{\ρ}A)=0---\left(3\right)$

求解上述方程得到ω_n和λ_n(λ_1n,λ_2n)之间的关系：

${\mathrm{\λ}}_{1,2n}=\sqrt{\frac{1}{2EI}\[{{\mathrm{\ω}}_n}^2\mathrm{\ρ}(\frac{EI}{\mathrm{\κ}G}+I)\±\sqrt{{{\mathrm{\ω}}_n}^4{\mathrm{\ρ}}^2{(\frac{EI}{\mathrm{\κ}G}+I)}^2-4EI(\frac{{{\mathrm{\ω}}_n}^4{\mathrm{\ρ}}^{2}I}{\mathrm{\κ}G}-{{\mathrm{\ω}}_n}^2\mathrm{\ρ}A)}\]}---\left(4\right)$

将ω_n和λ_n(λ_1n,λ_2n)之间的关系代入位移v_n(x)和转角φ_n(x)的标准模态函数得到：

v_n(x)＝c_1nch(λ_1nx)+s_1nsh(λ_1nx)+c_2ncos(λ_2nx)+s_2nsin(λ_2nx) (5)

φ_n(x)＝c_1ng_1nsh(λ_1nx)+s_1ng_1nch(λ_1nx)-c_2ng_2nsin(λ_2nx)+s_2ng_2ncos(λ_2nx) (6)

式中：

c_1n、c_2n、s_1n、s_2n为振型函数系数；

根据管节简化模型建立边界条件：

$\left\{\begin{array}{cc}\frac{d\mathrm{\φ}\left(x\right)}{dx}{|}_{x=0}=0,& \[\frac{dv\left(x\right)}{dx}-\mathrm{\φ}\left(x\right)\]{|}_{x=0}=0\\ \frac{d\mathrm{\φ}\left(x\right)}{dx}{|}_{x=l}=0,& \[\frac{dv\left(x\right)}{dx}-\mathrm{\φ}\left(x\right)\]{|}_{x=l}=0\end{array}\right.---\left(7\right)$

式中：l为管节长度；

满足模态函数系数c_1n、c_2n、s_1n、s_2n不同时等于0，求解管节振动固有频率ω_n，从而得到管节模态振型，具体采用Matlab编程求解；上述方法适用于弹性体模态求解，而根据相关研究^[17]，自由边界条件下Timoshenko梁前两阶模态为刚体模态，其模态函数及频率为：

$\left\{\begin{array}{ccc}{v}_1\left(x\right)=1,& {\mathrm{\φ}}_1\left(x\right)=0,& {\mathrm{\ω}}_1=0\\ v_2\left(x\right)=\sqrt{3}(1-2x/l),& {\mathrm{\φ}}_1\left(x\right)=0,& {\mathrm{\ω}}_1=0\end{array}\right.---\left(8\right)$

步骤2)：管节动力方程建立及求解

先建立管节受迫振动控制方程：

$\left\{\begin{array}{c}\mathrm{\κ}AG(\frac{{\∂}^{2}v(x,t)}{\∂x^2}-\frac{\∂\mathrm{\φ}(x,t)}{\∂x})=\mathrm{\ρ}A\frac{{\∂}^{2}v(x,t)}{\∂t^2}+F(x,t)\\ EI\frac{{\∂}^2\mathrm{\φ}(x,t)}{\∂x^2}+\mathrm{\κ}AG(\frac{\∂v(x,t)}{\∂x}-\mathrm{\φ}(x,t\left)\right)=\mathrm{\ρ}I\frac{{\∂}^2\mathrm{\φ}(x,t)}{\∂t^2}+M(x,t)\end{array}\right.---\left(9\right)$

式中：F(x,t)为管节所受外力，单位N/m；

M(x,t)为管节所受弯矩N·m/m；

采用模态叠加法，假定梁的竖向位移及转角表达式为：

$\left\{\begin{array}{c}v(x,t)=\underset{n=1}{\overset{me}{\mathrm{\Σ}}}{v}_n\left(x\right)q_n\left(t\right)\\ \mathrm{\φ}(x,t)=\underset{n=1}{\overset{me}{\mathrm{\Σ}}}{\mathrm{\φ}}_n\left(x\right)q_n\left(t\right)\end{array}\right.---\left(10\right)$

式中：q_n(t)为时间系数，单位为s；

将(10)代入(9)，进行正交化解耦得到第j段管节的第n阶振动常微分方程为：

${\∫}_0^{l_j}\left({\mathrm{\ρAv}}_n^{j2}\right(x)+{\mathrm{\ρI\φ}}_n^{j2}(x\left)\right)dx\·{\stackrel{\·\·}{q}}_n^j\left(t\right)+{\mathrm{\ω}}_n^{j2}{\∫}_0^{l_j}\left({\mathrm{\ρAv}}_n^{j2}\right(x)+{\mathrm{\ρI\φ}}_n^{j2}(x\left)\right)dx\·q_n^j\left(t\right)={\∫}_0^{l_j}{v}_n^j\left(x\right)F^j(x,t)dx+{\∫}_0^{l_j}{\mathrm{\φ}}_n^j\left(x\right)M^j(x,t)dx---\left(11\right)$

式中：l_j为第j段管节长度，单位为m；

由于车辆质量相对管节质量可忽略不计，本文将车辆前后轴荷载等效为两个点源移动恒载：

P(t)＝∑P_mδ(x-(ut+x_m))δ(y) (12)

式中：P_m为时刻第m辆车作用在管节上的点荷载，单位为N；

δ(·)为狄拉克函数；

u为车辆行驶速度，单位为m/s；

x_m为车辆初始位置，单位为m；

y为管节横向坐标，单位为m；

假设车辆沿隧道轴线方向行驶，考虑车辆荷载、地基反力和接头集中力及弯矩作用，Fⁱ(x,t)和Mⁱ(x,t)的具体表达式为：

$F^j\left(x,t\right)=\left\{\begin{array}{c}{\mathrm{\ΣP}}_{my}\mathrm{\δ}\left(x-\left(ut+x_m\right)\right)-{\mathrm{kv}}^j\left(x,t\right)-c{\stackrel{\·}{v}}^j\left(x,t\right)+\\ c_j\[\left({\stackrel{\·}{v}}^{j-1}\left(l_{j-1},t\right)-{\stackrel{\·}{v}}^j\left(0,t\right)\right)\mathrm{\δ}\left(x\right)+\left({\stackrel{\·}{v}}^{j+1}\left(0,t\right)-{\stackrel{\·}{v}}^j\left(l_j,t\right)\right)\mathrm{\δ}\left(x-l_j\right)\]+\\ k_j\[\left(v^{j-1}\left(l_{j-1},t\right)-v^j\left(0,t\right)\right)\mathrm{\δ}\left(x\right)+\left(v^{j+1}\left(0,t\right)-v^j\left(l_j,t\right)\right)\mathrm{\δ}\left(x-l_j\right)\]\end{array}\right\}---\left(13\right)$

式中：k_j为接头抗剪单元弹簧系数，单位为N/m；

c_j为接头抗剪单元阻尼系数N·s/m；

k为地基等效弹簧系数，N/m²；

c为地基阻尼系数N·s/m²；

P_my为车辆等效横向均布荷载，单位为N/m；

$M^j(x,t)=\left\{\begin{array}{c}{c}_w\[(\frac{{\∂}^2{v}^{j-1}}{\∂x\∂t}{|}_{x=l_{j-1}}-\frac{{\∂}^2{v}^j}{\∂x\∂t}{|}_{x=0})\mathrm{\δ}\left(x\right)+(\frac{{\∂}^2{v}^{j+1}}{\∂x\∂t}{|}_{x=0}-\frac{{\∂}^2{v}^j}{\∂x\∂t}{|}_{x=l_j})\mathrm{\δ}(x-l_j)\]+\\ k_w\[(\frac{\∂v^{j-1}}{\∂x}{|}_{x=l_{j-1}}-\frac{\∂v^j}{\∂x}{|}_{x=0})\mathrm{\δ}\left(x\right)+(\frac{\∂v^{j+1}}{\∂x}{|}_{x=0}-\frac{\∂v^j}{\∂x}{|}_{x=l_j})\mathrm{\δ}(x-l_j)\]\end{array}\right\}---\left(14\right)$

式中：k_w为接头抗弯单元弹簧系数，单位为N·m/rad；

c_w为接头抗弯单元阻尼系数，单位为N·m·s/rad；

将(13)和(14)代入(11)最终得到：

$\begin{array}{c}{\∫}_0^{l_j}\left({\mathrm{\ρAv}}_n^{j2}\left(x\right)+{\mathrm{\ρI\φ}}_n^{j2}\left(x\right)\right)dx\·{\stackrel{\·\·}{q}}_n^j\left(t\right)+{\mathrm{\ω}}_n^{j2}{\∫}_0^l\left({\mathrm{\ρAv}}_n^{j2}\left(x\right)+{\mathrm{\ρI\φ}}_n^{j2}\left(x\right)\right)dx\·q_n^j\left(t\right)\\ =\left\{\begin{array}{c}{\mathrm{\Σv}}_n^j\left(ut+x_m\right)\×P_{my}-\left({\mathrm{kq}}_n^j\left(t\right)+c{\stackrel{\·}{q}}_n^j\left(t\right)\right)\×{\∫}_0^{l_j}{v}_n^{j2}\left(x\right)dx+\\ v_n^j\left(0\right)\×\[c_j\×\left(\underset{u=1}{\overset{me}{\mathrm{\Σ}}}{v}_u^{j-1}\left(l_{j-1}\right){\stackrel{\·}{q}}_u^{j-1}\left(t\right)-\underset{u=1}{\overset{me}{\mathrm{\Σ}}}{v}_u^j\left(0\right){\stackrel{\·}{q}}_u^j\left(t\right)\right)+k_j\×\left(\underset{u=1}{\overset{me}{\mathrm{\Σ}}}{v}_u^{j-1}\left(l_{j-1}\right)q_u^{j-1}\left(t\right)-\underset{u=1}{\overset{me}{\mathrm{\Σ}}}{v}_u^j\left(0\right)q_u^j\left(t\right)\right)\]+\\ v_n^j\left(l_j\right)\×\[c_j\×\left(\underset{u=1}{\overset{me}{\mathrm{\Σ}}}{v}_u^{j+1}\left(0\right){\stackrel{\·}{q}}_u^{j+1}\left(t\right)-\underset{u=1}{\overset{me}{\mathrm{\Σ}}}{v}_u^j\left(l_j\right){\stackrel{\·}{q}}_u^j\left(t\right)\right)+k_j\×\left(\underset{u=1}{\overset{me}{\mathrm{\Σ}}}{v}_u^{j+1}\left(0\right)q_u^{j+1}\left(t\right)-\underset{u=1}{\overset{me}{\mathrm{\Σ}}}{v}_u^j\left(l_j\right)q_u^j\left(t\right)\right)\]+\\ {\mathrm{\φ}}_n^j\left(0\right)\×\[c_w\×\left(\underset{u=1}{\overset{me}{\mathrm{\Σ}}}{\stackrel{\·}{v}}_u^{j-1}\left(l_{j-1}\right){\stackrel{\·}{q}}_u^{j-1}\left(t\right)-\underset{u=1}{\overset{me}{\mathrm{\Σ}}}{\stackrel{\·}{v}}_u^j\left(0\right){\stackrel{\·}{q}}_u^j\left(t\right)\right)+k_w\×\left(\underset{u=1}{\overset{me}{\mathrm{\Σ}}}{\stackrel{\·}{v}}_u^{j-1}\left(l_{j-1}\right)q_u^{j-1}\left(t\right)-\underset{u=1}{\overset{me}{\mathrm{\Σ}}}{\stackrel{\·}{v}}_u^j\left(0\right)q_u^j\left(t\right)\right)\]+\\ {\mathrm{\φ}}_n^j\left(l_j\right)\×\[c_w\×\left(\underset{u=1}{\overset{me}{\mathrm{\Σ}}}{\stackrel{\·}{v}}_u^{j+1}\left(0\right){\stackrel{\·}{q}}_u^{j+1}\left(t\right)-\underset{u=1}{\overset{me}{\mathrm{\Σ}}}{\stackrel{\·}{v}}_u^j\left(l_j\right){\stackrel{\·}{q}}_u^j\left(t\right)\right)+k_w\×\left(\underset{u=1}{\overset{me}{\mathrm{\Σ}}}{\stackrel{\·}{v}}_u^{j+1}\left(0\right)q_u^{j+1}\left(t\right)-\underset{u=1}{\overset{me}{\mathrm{\Σ}}}{\stackrel{\·}{v}}_u^j\left(l_j\right)q_u^j\left(t\right)\right)\]\end{array}\right\}\end{array}---\left(15\right)$

将(15)整理成矩阵方程组，采用Newmark逐步积分法求解，得到第j段管节第n阶时间系数结合管节模态函数能够得到管节纵向任意位置的竖向位移响应。