强度计算.结构分析：动力学分析：有限元方法在动力学分析中的应用

强度计算.结构分析：动力学分析：有限元方法在动力学分析中的应用1绪论1.1有限元方法简介有限元方法（FiniteElementMethod,FEM）是一种数值分析技术，广泛应用于工程和科学领域，特别是结构分析、流体动力学、热传导和电磁学等。它将复杂的连续体结构分解为有限数量的简单单元，即“有限元”，通过在每个单元上应用数学模型，可以近似求解这些单元的物理行为，进而得到整个结构的解。这种方法能够处理形状复杂、边界条件多变的结构问题，是现代工程分析不可或缺的工具。1.1.1原理有限元方法基于变分原理和加权残值法。它首先将结构的连续域离散化为有限个单元，每个单元用一组节点来表示。然后，通过在每个节点上定义位移，将结构的连续问题转化为离散问题。接下来，利用单元的刚度矩阵和质量矩阵，建立整个结构的方程组，通过求解这些方程组，可以得到结构在不同载荷下的响应。1.1.2示例假设我们有一个简单的梁，需要使用有限元方法分析其在载荷作用下的变形。我们可以将梁离散化为多个线性单元，每个单元有两个节点，分别定义了节点的位移。下面是一个使用Python和SciPy库进行简单梁分析的代码示例：importnumpyasnp

fromscipy.sparseimportdiags

fromscipy.sparse.linalgimportspsolve

#定义梁的参数

E=200e9#弹性模量，单位：Pa

I=0.05**4/12#惯性矩，单位：m^4

L=1.0#梁的长度，单位：m

n_elements=10#单元数量

n_nodes=n_elements+1#节点数量

element_length=L/n_elements#单元长度

force=-1000#载荷，单位：N

#计算单元刚度矩阵

k=(E*I)/element_length**3*np.array([[12,6*element_length,-12,6*element_length],

[6*element_length,4*element_length**2,-6*element_length,2*element_length**2],

[-12,-6*element_length,12,-6*element_length],

[6*element_length,2*element_length**2,-6*element_length,4*element_length**2]])

#组装全局刚度矩阵

K=diags([np.repeat(k[0,0],n_elements),

np.repeat(k[0,1]+k[1,0],n_elements-1),

np.repeat(k[1,1],n_elements),

np.repeat(k[1,2]+k[2,1],n_elements-1),

np.repeat(k[2,2],n_elements),

np.repeat(k[2,3]+k[3,2],n_elements-1),

np.repeat(k[3,3],n_elements)],

[0,1,2,3,4,5,6],shape=(2*n_nodes,2*n_nodes)).toarray()

#定义边界条件和载荷向量

boundary_conditions=np.zeros(2*n_nodes)

boundary_conditions[0]=1#固定端

boundary_conditions[-1]=1#固定端

F=np.zeros(2*n_nodes)

F[n_nodes-1]=force#应用载荷

#应用边界条件

K[0,:]=0

K[-1,:]=0

K[:,0]=0

K[:,-1]=0

K[0,0]=1

K[-1,-1]=1

#求解位移向量

U=spsolve(K,F)

#输出结果

print("节点位移：",U)1.2动力学分析的重要性动力学分析是结构分析的一个重要分支，它研究结构在动态载荷作用下的响应，如振动、冲击和地震等。动力学分析对于确保结构在各种动态条件下的安全性和稳定性至关重要。例如，桥梁、建筑物和飞机等结构在设计时必须考虑风、地震和飞行过程中的动态载荷，以避免结构的破坏或失效。1.2.1应用动力学分析在多个领域有广泛的应用：-航空航天：分析飞机在飞行过程中的振动和稳定性。-土木工程：评估桥梁和建筑物在地震、风和交通载荷下的响应。-机械工程：研究机器部件在运行过程中的振动和疲劳。1.3结构动力学基础理论结构动力学主要研究结构在时间变化载荷作用下的响应。它涉及到牛顿第二定律的应用，即力等于质量乘以加速度。在有限元方法中，结构的动力学响应可以通过求解运动方程来获得，运动方程通常表示为：M其中，M是质量矩阵，C是阻尼矩阵，K是刚度矩阵，u是位移向量，u和u分别表示位移的一阶和二阶导数，Ft1.3.1解析方法解析方法通常用于解决线性动力学问题，其中结构的响应可以通过求解微分方程的解析解来获得。然而，对于复杂结构和非线性问题，解析方法往往难以应用，这时就需要使用数值方法，如有限元方法。1.3.2数值方法数值方法，如有限元方法，通过将结构离散化为有限个单元，然后在每个单元上应用运动方程，可以求解结构的动力学响应。这种方法能够处理非线性问题和复杂边界条件，是现代工程分析中解决动力学问题的主要手段。1.3.3示例下面是一个使用Python和SciPy库进行结构动力学分析的代码示例，分析一个简单的单自由度系统在冲击载荷作用下的响应：importnumpyasnp

fromegrateimportsolve_ivp

#定义系统参数

m=1.0#质量，单位：kg

k=1000#弹性系数，单位：N/m

c=10#阻尼系数，单位：N/(m/s)

#定义冲击载荷

defforce(t):

ift<0.1:

return10000#冲击载荷，单位：N

else:

return0

#定义运动方程

defmotion(t,y):

u,v=y#位移和速度

du_dt=v

dv_dt=(-c*v-k*u+force(t))/m

return[du_dt,dv_dt]

#初始条件

y0=[0,0]#初始位移和速度

#时间范围

t_span=(0,1)

#求解运动方程

sol=solve_ivp(motion,t_span,y0,t_eval=np.linspace(0,1,100))

#输出结果

importmatplotlib.pyplotasplt

plt.plot(sol.t,sol.y[0],label='Displacement')

plt.plot(sol.t,sol.y[1],label='Velocity')

plt.legend()

plt.show()这个例子展示了如何使用数值积分方法（如solve_ivp）来求解单自由度系统的运动方程，得到位移和速度随时间变化的曲线。2有限元方法基础2.1节点与单元的概念在有限元分析中，结构被离散化为一系列小的、可管理的单元，这些单元通过节点连接在一起。节点是结构中的特定点，它们是单元的边界，也是力和位移的计算点。单元则是结构的组成部分，可以是线、面或体，它们的形状和大小取决于分析的精度和计算资源。2.1.1节点节点是有限元模型中的基本元素，它们定义了结构的几何形状和边界条件。在每个节点上，可以定义位移、速度、加速度等自由度，以及作用在节点上的外力。节点之间的连接关系决定了单元的形状和类型。2.1.2单元单元是连接节点的几何体，它们可以是线性的（如梁、杆）、平面的（如壳、板）或三维的（如实体）。单元的类型决定了其内部应力和应变的计算方式。例如，实体单元通常用于分析三维结构的应力分布，而壳单元则用于薄壁结构的分析。2.2材料属性与边界条件材料属性和边界条件是有限元分析中不可或缺的部分，它们直接影响分析结果的准确性和可靠性。2.2.1材料属性材料属性包括弹性模量、泊松比、密度等，这些属性用于描述材料在受力时的行为。例如，弹性模量反映了材料抵抗弹性变形的能力，泊松比描述了材料在拉伸或压缩时横向变形的特性。2.2.2边界条件边界条件定义了结构的约束和外部作用。常见的边界条件包括固定约束、滑动约束、力和力矩的施加等。正确设置边界条件对于模拟真实结构的响应至关重要。2.3有限元模型的建立建立有限元模型是一个系统的过程，包括几何建模、网格划分、材料属性和边界条件的定义、以及模型的求解和后处理。2.3.1几何建模首先，需要在CAD软件中创建结构的几何模型。这一步骤确保了模型的几何精度，是后续分析的基础。2.3.2网格划分网格划分是将几何模型离散化为单元和节点的过程。网格的大小和形状需要根据分析的精度要求和计算资源来确定。细网格可以提供更精确的分析结果，但会增加计算时间和资源需求。2.3.3材料属性和边界条件的定义在模型中定义每个单元的材料属性和所有节点的边界条件。这一步骤需要根据实际结构的材料和工作环境来设定。2.3.4模型求解使用有限元软件对模型进行求解，计算在给定载荷下的结构响应，如位移、应力和应变。2.3.5后处理分析求解结果，可视化应力分布、位移等，以评估结构的性能和安全性。2.3.6示例：使用Python和FEniCS进行有限元分析fromdolfinimport*

#创建一个矩形网格

mesh=RectangleMesh(Point(0,0),Point(1,1),10,10)

#定义边界条件

defboundary(x,on_boundary):

returnon_boundary

bc=DirichletBC(VectorFunctionSpace(mesh,'CG',1),Constant((0,0)),boundary)

#定义材料属性

E=1e3#弹性模量

nu=0.3#泊松比

mu=E/(2*(1+nu))

lmbda=E*nu/((1+nu)*(1-2*nu))

#定义变分问题

V=VectorFunctionSpace(mesh,'CG',1)

u=TrialFunction(V)

v=TestFunction(V)

f=Constant((0,-10))#作用力

T=Constant(0)#温度

#应力应变关系

defsigma(u):

returnlmbda*tr(eps(u))*Identity(2)+2.0*mu*eps(u)

#应变位移关系

defeps(u):

returnsym(nabla_grad(u))

#定义变分形式

a=inner(sigma(u),eps(v))*dx

L=dot(f,v)*dx+dot(T,v)*ds

#求解

u=Function(V)

solve(a==L,u,bc)

#后处理

plot(u)

interactive()在这个示例中，我们使用Python的FEniCS库来建立和求解一个简单的有限元模型。模型是一个1x1的矩形，受到垂直向下的力作用。我们定义了边界条件，材料属性（弹性模量和泊松比），并使用变分原理来建立和求解力学问题。最后，我们通过可视化位移来检查分析结果。3动力学分析类型在结构分析领域，动力学分析是评估结构在动态载荷作用下行为的关键步骤。本教程将深入探讨三种主要的动力学分析类型：瞬态动力学分析、模态分析和谐波响应分析，以及有限元方法在这些分析中的应用。3.1瞬态动力学分析3.1.1原理瞬态动力学分析关注结构在时间域内的响应，特别是当结构受到随时间变化的载荷作用时。这种分析可以预测结构在特定时间点的位移、速度、加速度和应力，对于理解结构在地震、爆炸、冲击等瞬态事件下的行为至关重要。3.1.2内容瞬态动力学分析通常基于牛顿第二定律，即力等于质量乘以加速度。在有限元分析中，结构被离散成多个小的单元，每个单元的运动方程被建立并求解。这些方程通常是非线性的，需要使用数值方法，如Newmark方法或中央差分法，来求解。示例假设我们有一个简单的弹簧-质量系统，质量为1kg，弹簧刚度为100N/m，受到一个随时间变化的力的作用。我们可以使用Python的egrate.solve_ivp函数来求解系统的运动方程。importnumpyasnp

fromegrateimportsolve_ivp

#定义系统运动方程

defspring_mass(t,y,k,m):

dydt=[y[1],-k/m*y[0]]

returndydt

#参数

k=100#弹簧刚度

m=1#质量

y0=[0,0]#初始条件：位移和速度

#时间范围

t_span=(0,10)

#求解

sol=solve_ivp(spring_mass,t_span,y0,args=(k,m),t_eval=np.linspace(0,10,100))

#输出结果

print(sol.t)#时间点

print(sol.y[0])#位移3.2模态分析3.2.1原理模态分析用于确定结构的固有频率、振型和阻尼比。这些信息对于理解结构的动态特性至关重要，特别是在设计阶段，以避免共振和优化结构的动态性能。3.2.2内容模态分析基于结构的自由振动方程，即在没有外部载荷作用下，结构的振动方程。在有限元分析中，结构的刚度矩阵和质量矩阵被用来求解这些方程，从而得到固有频率和振型。示例使用Python的scipy.linalg.eig函数，我们可以求解一个结构的固有频率和振型。假设我们有一个由两个质量块和三个弹簧组成的系统。importnumpyasnp

fromscipy.linalgimporteig

#定义质量矩阵和刚度矩阵

M=np.array([[1,0],[0,1]])#质量矩阵

K=np.array([[200,-100],[-100,200]])#刚度矩阵

#求解固有频率和振型

eigenvalues,eigenvectors=eig(-np.linalg.inv(M)*K)

#固有频率为固有值的平方根

frequencies=np.sqrt(eigenvalues)

#输出结果

print(frequencies)#固有频率

print(eigenvectors)#振型3.3谐波响应分析3.3.1原理谐波响应分析用于评估结构在正弦载荷作用下的响应。这种分析特别适用于预测结构在周期性载荷下的行为，如旋转机械的振动。3.3.2内容在谐波响应分析中，结构的响应被假设为正弦函数，其频率与载荷的频率相同。通过求解结构的运动方程，可以得到结构在不同频率下的响应，包括位移、速度和加速度。示例使用Python的scipy.signal库，我们可以模拟一个结构在谐波载荷作用下的响应。假设我们有一个质量为1kg，弹簧刚度为100N/m的系统，受到频率为1Hz的正弦载荷作用。importnumpyasnp

fromscipy.signalimportlti,impulse

#定义系统参数

k=100#弹簧刚度

m=1#质量

system=lti([1],[m,0,k])#创建一个连续时间LTI系统

#定义输入信号：频率为1Hz的正弦波

t=np.linspace(0,10,1000)

u=np.sin(2*np.pi*t)

#求解系统响应

t,y,x=impulse(system,T=t)

#输出结果

print(t)#时间点

print(y)#位移响应通过以上示例，我们可以看到，有限元方法在动力学分析中的应用不仅限于理论计算，还可以通过编程实现具体问题的求解，为工程设计和分析提供强大的工具。4有限元软件操作4.1前处理：模型建立与网格划分在进行动力学分析之前，前处理阶段是至关重要的，它包括模型的建立和网格的划分。这一阶段决定了分析的准确性和效率。4.1.1模型建立模型建立涉及定义结构的几何形状、材料属性和边界条件。例如，假设我们正在分析一个简单的梁结构，其材料为钢，弹性模量为200GPa，泊松比为0.3。4.1.2网格划分网格划分是将连续的结构离散化为有限数量的单元。单元的大小和形状对分析结果有直接影响。对于复杂的结构，可能需要使用非结构化网格。#以下是一个使用Python和FEniCS库进行网格划分的示例

fromfenicsimport*

#创建一个矩形网格

mesh=RectangleMesh(Point(0,0),Point(1,0.1),10,1)

#定义函数空间

V=FunctionSpace(mesh,'P',1)

#定义边界条件

defboundary(x,on_boundary):

returnon_boundary

bc=DirichletBC(V,Constant(0),boundary)

#定义材料属性

E=200e9#弹性模量

nu=0.3#泊松比

rho=7800#密度

#定义方程

#...（此处省略具体方程定义）4.2求解设置：加载与求解器选择在动力学分析中，加载和求解器的选择是关键步骤。加载可以是静态的，也可以是动态的，如冲击载荷或振动。4.2.1加载加载定义了作用在结构上的力或位移。例如，我们可以定义一个冲击载荷，其在短时间内施加一个大的力。4.2.2求解器选择求解器的选择取决于问题的性质。对于线性问题，可以直接使用直接求解器。对于非线性问题，可能需要迭代求解器。#以下是一个使用Python和FEniCS库进行加载和求解器选择的示例

fromfenicsimport*

#定义加载

f=Constant((0,-10))

#定义方程

#...（此处省略具体方程定义）

#定义求解器

solver=LinearSolver('mumps')

#求解问题

#...（此处省略具体求解步骤）4.3后处理：结果可视化与分析后处理阶段涉及分析和可视化求解结果。这可以帮助我们理解结构的动力学行为，如应力分布、位移和振动模式。4.3.1结果可视化结果可视化通常使用专门的软件，如Paraview或Mayavi。然而，FEniCS库也提供了基本的可视化功能。4.3.2结果分析结果分析可能包括计算结构的模态、频率响应或瞬态响应。例如，我们可以计算梁的前几阶模态。#以下是一个使用Python和FEniCS库进行结果可视化和分析的示例

fromfenicsimport*

importmatplotlib.pyplotasplt

#定义结果

u=Function(V)

#读取结果

#...（此处省略读取结果的步骤）

#可视化结果

plot(u)

plt.show()

#分析结果

#...（此处省略具体分析步骤，如计算模态）通过以上步骤，我们可以使用有限元方法进行动力学分析，从模型建立到网格划分，再到求解设置和后处理，每一步都对最终的分析结果至关重要。5动力学分析案例研究5.1桥梁结构的动力学分析5.1.1原理与内容桥梁结构的动力学分析是评估桥梁在动态载荷作用下（如风、地震、车辆通过等）的响应和稳定性。有限元方法（FEM）是解决此类问题的常用工具，它将桥梁结构离散成多个小的、简单的单元，每个单元的力学行为可以被精确描述，然后通过组合这些单元来模拟整个结构的动态响应。有限元模型建立几何建模：使用CAD软件创建桥梁的三维模型。材料属性：定义桥梁各部分的材料属性，如混凝土、钢材的弹性模量、泊松比等。网格划分：将模型划分为足够小的单元，以确保分析的准确性。动态载荷施加地震载荷：通过输入地震波或使用模态分析来模拟地震效应。风载荷：考虑风速和风向，使用流体动力学分析来模拟风对桥梁的影响。车辆载荷：模拟车辆通过桥梁时的动态效应，包括车辆重量、速度和分布。分析与结果模态分析：确定桥梁的固有频率和振型。瞬态分析：计算桥梁在特定动态载荷下的时间响应。谐响应分析：评估桥梁在周期性载荷下的响应。5.1.2示例：桥梁模态分析假设我们有一个简化的桥梁模型，使用Python的scipy库进行模态分析。importnumpyasnp

fromscipy.linalgimporteig

#定义桥梁的刚度矩阵和质量矩阵

K=np.array([[4,-2],[-2,4]])#刚度矩阵

M=np.array([[2,0],[0,2]])#质量矩阵

#计算固有频率和振型

eigenvalues,eigenvectors=eig(-np.linalg.inv(M)@K)

#转换为Hz

frequencies=np.sqrt(eigenvalues)/(2*np.pi)

#输出结果

print("固有频率：",frequencies)

print("振型：",eigenvectors)此代码示例中，我们使用了两个自由度的简化桥梁模型，通过计算刚度矩阵和质量矩阵的特征值和特征向量来确定固有频率和振型。实际应用中，桥梁模型将更加复杂，包含成千上万个自由度。5.2高层建筑的地震响应分析5.2.1原理与内容高层建筑的地震响应分析旨在评估建筑在地震作用下的安全性和性能。有限元方法通过模拟地震波对建筑的影响，计算结构的位移、加速度和内力，以确保建筑能够承受预期的地震载荷。有限元模型建立结构建模：创建建筑的三维模型，包括基础、柱、梁和楼板。材料属性：定义混凝土、钢材等的力学性能。边界条件：设置基础的固定约束，模拟地基的刚性。地震载荷施加地震波输入：使用实际记录的地震波或合成地震波作为输入。地震效应：考虑地震波的频率、振幅和方向，模拟地震对建筑的影响。分析与结果非线性动力分析：考虑材料的非线性行为和结构的大变形。损伤评估：评估地震后结构的损伤程度。安全评估：确保建筑满足抗震设计规范。5.2.2示例：地震响应分析使用Python的PyDSTool库进行地震响应分析的简化示例。fromPyDSToolimport*

fromPyDSTool.Toolboximportphaser

#定义建筑模型参数

params={'m':1000,'k':100000,'c':100}

#创建模型

model=ODEsystem('model',{'m':'m','k':'k','c':'c'},

{'x':'x','v':'v'},

{'dxdt':'v','dvdt':'-k*x/m-c*v/m'},

params)

#设置初始条件和时间范围

ICs={'x':0,'v':0}

t=np.linspace(0,10,1000)

#输入地震波

earthquake=np.sin(2*np.pi*1*t)#简化地震波

#进行动态分析

traj=pute('traj',ICs,t,{'earthquake':earthquake})

#可视化结果

phaser.plotSolution(traj,'x','v')此示例中，我们创建了一个单自由度的建筑模型，并施加了一个简化的正弦波作为地震载荷。实际分析中，将使用更复杂的模型和真实的地震波数据。5.3旋转机械的振动分析5.3.1原理与内容旋转机械的振动分析用于识别和预测旋转部件（如轴、齿轮、涡轮机等）在运行过程中的振动特性。有限元方法可以模拟旋转机械的动态行为，帮助工程师优化设计，减少振动和噪音。有限元模型建立几何建模：创建旋转机械的三维模型。材料属性：定义材料的弹性模量、密度等。旋转条件：设置旋转速度和方向。动态载荷施加不平衡载荷：模拟旋转部件的不平衡效应。摩擦载荷：考虑部件之间的摩擦力。流体载荷：如果适用，模拟流体对旋转部件的影响。分析与结果频谱分析：识别振动的主要频率。振幅分析：计算振动的振幅。模式分析：确定振动模式。5.3.2示例：旋转轴的不平衡振动分析使用Python的numpy和matplotlib库进行不平衡振动分析的简化示例。importnumpyasnp

importmatplotlib.pyplotasplt

#定义旋转轴的参数

mass=1.0#质量

radius=0.1#半径

unbalance=0.01#不平衡量

omega=2*np.pi*100#角速度

#计算振动响应

time=np.linspace(0,1,1000)

displacement=unbalance*radius*np.cos(omega*time)

#绘制振动响应

plt.figure()

plt.plot(time,displacement)

plt.title('旋转轴的不平衡振动响应')

plt.xlabel('时间(s)')

plt.ylabel('位移(m)')

plt.show()此代码示例中，我们模拟了一个旋转轴在不平衡载荷作用下的振动响应。实际应用中，将考虑更复杂的旋转机械系统，并使用有限元软件进行详细的振动分析。以上案例研究展示了有限元方法在动力学分析中的应用，从桥梁结构到高层建筑，再到旋转机械，每种结构都有其特定的分析需求和挑战。通过建立准确的有限元模型，施加适当的动态载荷，并进行详细的分析，工程师可以确保结构的安全性和性能。6高级主题与挑战6.1非线性动力学分析6.1.1原理非线性动力学分析考虑了结构在大变形、大应变、材料非线性、接触非线性等条件下的行为。在有限元方法中，非线性问题的求解通常需要迭代过程，以逐步逼近真实解。非线性动力学分析的关键在于正确处理非线性项，确保在每个时间步长内求解的稳定性与收敛性。6.1.2内容材料非线性：如塑性、粘弹性、超弹性等。几何非线性：大变形、大应变效应。接触非线性：结构间接触、摩擦等。边界条件非线性：非线性弹簧、阻尼器等。6.1.3示例：材料非线性分析假设有一个简单的弹簧-质量系统，其中弹簧表现出非线性行为，其力-位移关系由以下方程描述：F其中，F是力，u是位移，k是非线性刚度系数。Python代码示例importnumpyasnp

fromegrateimportsolve_ivp

#定义非线性弹簧的力-位移关系

defnonlinear_force(u,k):

returnk*u**3

#定义动力学方程

defdynamics(t,y,k,m):

u,v=y

du_dt=v

dv_dt=-nonlinear_force(u,k)/m

return[du_dt,dv_dt]

#参数设置

m=1.0#质量

k=1.0#初始刚度系数

t_span=(0,10)#时间跨度

y0=[0.1,0.0]#初始条件：位移和速度

#求解动力学方程

sol=solve_ivp(dynamics,t_span,y0,args=(k,m),t_eval=np.linspace(t_span[0],t_span[1],100))

#输出结果

importmatplotlib.pyplotasplt

plt.plot(sol.t,sol.y[0],label='位移')

plt.plot(sol.t,sol.y[1],label='速度')

plt.legend()

plt.show()解释此代码示例使用Python的egrate.solve_ivp函数求解非线性动力学方程。dynamics函数定义了系统的动力学方程，其中考虑了非线性弹簧的力-位移关系。通过调整k值，可以模拟不同非线性程度的弹簧。solve_ivp函数用于数值积分，求解系统的位移和速度随时间的变化。6.2多物理场耦合动力学分析6.2.1原理多物理场耦合动力学分析涉及结构动力学与热、电、磁等其他物理场的相互作用。在有限元分析中，这意味着需要同时求解多个物理场的方程，并考虑它们之间的耦合效应。例如，热效应可能改变材料的力学性能，从而影响结构的动力学响应。6.2.2内容热-结构耦合：温度变化引起材料性能变化，影响结构动力学。电-结构耦合：电磁力对结构的影响。流体-结构耦合：流体动力学对结构动力学的影响。6.2.3示例：热-结构耦合分析考虑一个受热的梁，其温度变化导致材料膨胀，从而影响梁的位移。Python代码示例importnumpyasnp

fromegrateimportsolve_ivp

#定义热-结构耦合方程

defcoupled_thermo_mechanics(t,y,k,m,alpha,q):

u,v,T=y

du_dt=v

dv_dt=-k*u/m+alpha*T

dT_dt=q/(m*1000)#假设比热为1000J/kg.K

return[du_dt,dv_dt,dT_dt]

#参数设置

m=1.0#质量

k=1.0#刚度

alpha=0.00