結(jié)構(gòu)力學數(shù)值方法：有限元法(FEM)：桿單元與梁單元分析

結(jié)構(gòu)力學數(shù)值方法：有限元法(FEM)：桿單元與梁單元分析1緒論1.1有限元法的基本概念有限元法（FiniteElementMethod,FEM）是一種用于求解復雜工程問題的數(shù)值分析方法。它將連續(xù)的結(jié)構(gòu)或系統(tǒng)離散成有限個單元，每個單元用簡單的數(shù)學模型來近似描述，然后通過組合這些單元的模型來求解整個系統(tǒng)的響應(yīng)。這種方法特別適用于解決那些無法通過解析方法求解的復雜結(jié)構(gòu)問題。1.1.1桿單元桿單元是最基本的有限元單元之一，主要用于分析一維結(jié)構(gòu)，如橋梁的拉桿或壓桿。桿單元的數(shù)學模型基于一維彈性方程，考慮軸向力的影響。在計算中，桿單元的長度、截面積、材料屬性（如彈性模量）以及兩端的約束條件是關(guān)鍵參數(shù)。1.1.2梁單元梁單元用于分析二維或三維結(jié)構(gòu)中的彎曲和剪切效應(yīng)，如橋梁的主梁或建筑的橫梁。梁單元的數(shù)學模型基于歐拉-伯努利梁理論或蒂蒙斯-維納梁理論，考慮彎矩、剪力和軸向力的影響。在計算中，梁單元的長度、截面積、慣性矩、材料屬性以及兩端的約束條件是關(guān)鍵參數(shù)。1.2結(jié)構(gòu)力學中的數(shù)值方法簡介結(jié)構(gòu)力學中的數(shù)值方法是解決結(jié)構(gòu)分析問題的一種有效手段，尤其當結(jié)構(gòu)形狀復雜、載荷分布不均或材料性質(zhì)非線性時，傳統(tǒng)的解析方法往往難以應(yīng)用。數(shù)值方法通過將結(jié)構(gòu)離散化，將其轉(zhuǎn)換為一系列可以計算的數(shù)學問題，從而能夠求解結(jié)構(gòu)的應(yīng)力、應(yīng)變和位移等關(guān)鍵參數(shù)。1.2.1離散化過程離散化是將連續(xù)的結(jié)構(gòu)轉(zhuǎn)換為有限個單元的過程。每個單元用節(jié)點來表示，節(jié)點之間的連接定義了單元的形狀和尺寸。通過在每個節(jié)點上應(yīng)用平衡方程，可以建立整個結(jié)構(gòu)的平衡方程組，進而求解結(jié)構(gòu)的響應(yīng)。1.2.2平衡方程組平衡方程組是有限元分析的核心。它基于牛頓第二定律，即力等于質(zhì)量乘以加速度。在結(jié)構(gòu)力學中，平衡方程組描述了結(jié)構(gòu)在外部載荷作用下的內(nèi)部力和位移之間的關(guān)系。通過求解這些方程，可以得到結(jié)構(gòu)在載荷作用下的響應(yīng)。1.3桿單元與梁單元的應(yīng)用領(lǐng)域桿單元和梁單元在結(jié)構(gòu)工程中有著廣泛的應(yīng)用，它們可以用于分析各種類型的結(jié)構(gòu)，包括但不限于：橋梁工程：分析橋梁的拉桿和主梁，評估其在不同載荷條件下的性能。建筑結(jié)構(gòu)：分析建筑中的柱、梁和桁架，確保結(jié)構(gòu)的安全性和穩(wěn)定性。機械工程：分析機械部件的應(yīng)力和應(yīng)變，如軸、連桿和梁，以優(yōu)化設(shè)計和提高可靠性。航空航天工程：分析飛機和航天器的結(jié)構(gòu)，確保其在極端條件下的安全性和性能。1.3.1示例：桿單元分析假設(shè)我們有一個簡單的桿單元，長度為1米，截面積為0.01平方米，材料的彈性模量為200GPa。兩端分別受到10kN的拉力和壓力。我們可以使用以下步驟來分析這個桿單元：定義單元參數(shù)：長度L=1m，截面積A建立平衡方程：對于桿單元，平衡方程可以簡化為σ=FA，其中σ求解應(yīng)力：將參數(shù)代入方程，得到σ=1.3.2示例：梁單元分析考慮一個簡單的梁單元，長度為2米，截面積為0.02平方米，慣性矩為0.0004m定義單元參數(shù)：長度L=2m，截面積A=0.02建立平衡方程：對于梁單元，平衡方程基于歐拉-伯努利梁理論，可以表示為d2dx2E求解撓度：在給定的邊界條件下，通過數(shù)值方法求解上述微分方程，可以得到梁在載荷作用下的撓度分布。通過這些基本概念和示例，我們可以看到桿單元和梁單元在結(jié)構(gòu)力學數(shù)值分析中的重要性和應(yīng)用價值。有限元法不僅能夠提供結(jié)構(gòu)的靜態(tài)分析，還可以進行動態(tài)分析、熱分析等多種類型的工程分析，是現(xiàn)代工程設(shè)計和分析不可或缺的工具。2有限元法基礎(chǔ)2.1離散化過程詳解在結(jié)構(gòu)力學的數(shù)值分析中，有限元法（FEM）是一種將連續(xù)體離散化為有限數(shù)量的單元進行分析的方法。這一過程通常包括以下步驟：幾何離散化：將結(jié)構(gòu)的幾何形狀分解為多個小的、簡單的幾何形狀，如桿、梁、板或?qū)嶓w單元。每個小單元可以視為結(jié)構(gòu)的一部分，其行為可以通過簡單的數(shù)學模型來描述。選擇位移函數(shù)：在每個單元內(nèi)，位移被假設(shè)為節(jié)點位移的函數(shù)。這些函數(shù)通常為多項式，能夠近似單元內(nèi)的實際位移分布。構(gòu)建單元方程：基于彈性力學的基本原理，如胡克定律和虛功原理，為每個單元建立方程。這些方程描述了單元內(nèi)部力與位移之間的關(guān)系。組裝整體方程：將所有單元方程組合成一個整體的結(jié)構(gòu)方程，通常表示為剛度矩陣方程。這個方程描述了整個結(jié)構(gòu)的力與位移之間的關(guān)系。施加邊界條件：在整體方程中加入結(jié)構(gòu)的邊界條件，如固定端、自由端或施加的外力和力矩。求解方程：使用數(shù)值方法，如直接求解或迭代求解，來求解整體方程，得到結(jié)構(gòu)的位移、應(yīng)力和應(yīng)變。2.1.1示例：桿單元的離散化與分析假設(shè)我們有一個簡單的桿結(jié)構(gòu)，長度為1米，兩端固定，中間受到1000N的集中力作用。我們將桿離散化為兩個單元，每個單元長度為0.5米。#Python示例代碼

#導入必要的庫

importnumpyasnp

#桿的屬性

E=200e9#彈性模量，單位：Pa

A=0.001#截面積，單位：m^2

L=1.0#總長度，單位：m

F=1000#集中力，單位：N

#單元長度

L1=L/2

L2=L/2

#單元剛度矩陣

K1=(E*A/L1)*np.array([[1,-1],[-1,1]])

K2=(E*A/L2)*np.array([[1,-1],[-1,1]])

#整體剛度矩陣

K=np.zeros((4,4))

K[0:2,0:2]+=K1

K[1:3,1:3]+=K2

K[2:4,2:4]+=K1

#邊界條件

#假設(shè)兩端固定，中間受力

boundary_conditions=np.array([1,0,1,0])

force=np.array([0,F,0,0])

#求解位移

displacements=np.linalg.solve(K[boundary_conditions==0,:][:,boundary_conditions==0],force[boundary_conditions==0])在這個例子中，我們首先定義了桿的物理屬性，包括彈性模量、截面積、總長度和作用力。然后，我們將桿離散化為兩個單元，并為每個單元構(gòu)建剛度矩陣。最后，我們組裝整體剛度矩陣，施加邊界條件，并求解位移。2.2節(jié)點與單元的概念在有限元分析中，節(jié)點是結(jié)構(gòu)中定義位移和力的點，而單元是連接這些節(jié)點的幾何體，用于近似結(jié)構(gòu)的物理行為。節(jié)點和單元是有限元模型的基本組成部分。節(jié)點：每個節(jié)點都有一個或多個自由度（DOF），如位移或旋轉(zhuǎn)。在桿單元中，每個節(jié)點通常有兩個自由度：橫向位移和軸向位移。單元：單元可以是線性的（如桿單元）或非線性的（如梁單元）。每個單元都有自己的剛度矩陣，描述了單元內(nèi)部力與位移之間的關(guān)系。2.2.1示例：梁單元的節(jié)點與單元考慮一個簡單的梁單元，長度為1米，兩端自由，中間受到1000N的集中力作用。我們將梁離散化為兩個單元，每個單元長度為0.5米。#Python示例代碼

#導入必要的庫

importnumpyasnp

#梁的屬性

E=200e9#彈性模量，單位：Pa

I=0.0001#慣性矩，單位：m^4

L=1.0#總長度，單位：m

F=1000#集中力，單位：N

#單元長度

L1=L/2

L2=L/2

#單元剛度矩陣

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

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

[-12,-6*L1,12,-6*L1],

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

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

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

[-12,6*L2,12,-6*L2],

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

#整體剛度矩陣

K=np.zeros((8,8))

K[0:4,0:4]+=K1

K[2:6,2:6]+=K2

K[4:8,4:8]+=K1

#邊界條件

#假設(shè)兩端自由，中間受力

boundary_conditions=np.array([0,0,1,0,0,0,1,0])

force=np.array([0,0,0,F,0,0,0,0])

#求解位移

displacements=np.linalg.solve(K[boundary_conditions==0,:][:,boundary_conditions==0],force[boundary_conditions==0])在這個例子中，我們定義了梁的物理屬性，包括彈性模量、慣性矩、總長度和作用力。然后，我們將梁離散化為兩個單元，并為每個單元構(gòu)建剛度矩陣。最后，我們組裝整體剛度矩陣，施加邊界條件，并求解位移。2.3剛度矩陣的構(gòu)建剛度矩陣是有限元分析的核心，它描述了結(jié)構(gòu)的力與位移之間的關(guān)系。對于不同的單元類型，剛度矩陣的構(gòu)建方法也不同。2.3.1桿單元的剛度矩陣對于桿單元，剛度矩陣通常為2x2矩陣，表示軸向力與軸向位移之間的關(guān)系。其形式為：K其中，E是彈性模量，A是截面積，L是單元長度。2.3.2梁單元的剛度矩陣對于梁單元，剛度矩陣通常為4x4矩陣，表示彎矩、剪力與轉(zhuǎn)角、橫向位移之間的關(guān)系。其形式為：K其中，E是彈性模量，I是慣性矩，L是單元長度。通過構(gòu)建和組裝這些剛度矩陣，我們可以得到整個結(jié)構(gòu)的力與位移之間的關(guān)系，從而進行結(jié)構(gòu)的數(shù)值分析。3桿單元分析3.1維桿單元的數(shù)學模型在結(jié)構(gòu)力學中，一維桿單元通常用于模擬承受軸向力的結(jié)構(gòu)，如橋梁的拉桿或建筑物的柱子。桿單元的數(shù)學模型基于胡克定律，即應(yīng)力與應(yīng)變成正比。對于一個簡單的桿單元，其長度為L，截面積為A，彈性模量為E，兩端分別受到軸向力F1和F3.1.1胡克定律σ其中，σ是應(yīng)力，?是應(yīng)變。3.1.2平衡方程F這里，u是位移，x是桿單元的坐標。3.1.3位移與力的關(guān)系F其中，k=3.2桿單元的剛度矩陣推導3.2.1剛度矩陣的定義剛度矩陣描述了結(jié)構(gòu)中力與位移之間的關(guān)系。對于一維桿單元，其剛度矩陣K是一個2×2的矩陣，表示兩端力F1和F2與兩端位移3.2.2推導過程考慮一個桿單元，兩端位移分別為u1和u2，則桿單元的總位移為u=u2?uF將兩端力表示為F1和FFF將上述方程寫成矩陣形式，我們得到桿單元的剛度矩陣：F3.2.3代碼示例#定義桿單元的剛度矩陣

defstiffness_matrix(A,E,L):

"""

計算一維桿單元的剛度矩陣

:paramA:截面積

:paramE:彈性模量

:paramL:桿單元長度

:return:剛度矩陣

"""

k=A*E/L

K=[[-k,k],[k,-k]]

returnK

#示例參數(shù)

A=0.01#截面積，單位：平方米

E=200e9#彈性模量，單位：帕斯卡

L=1.0#桿單元長度，單位：米

#計算剛度矩陣

K=stiffness_matrix(A,E,L)

print(K)3.3邊界條件的處理在有限元分析中，邊界條件的處理至關(guān)重要，它決定了結(jié)構(gòu)的約束方式。邊界條件可以分為位移邊界條件和力邊界條件。3.3.1位移邊界條件位移邊界條件通常用于固定結(jié)構(gòu)的某一部分，例如，將一端固定為u13.3.2力邊界條件力邊界條件用于施加外部力，例如，在一端施加力F13.3.3處理方法在計算過程中，邊界條件需要被應(yīng)用到剛度矩陣和力向量中。對于位移邊界條件，可以通過修改剛度矩陣和力向量來實現(xiàn)；對于力邊界條件，可以直接將力值添加到力向量中。3.3.4代碼示例#定義力向量

defforce_vector(F1,F2):

"""

計算一維桿單元的力向量

:paramF1:一端的力

:paramF2:另一端的力

:return:力向量

"""

F=[F1,F2]

returnF

#示例參數(shù)

F1=1000#一端的力，單位：牛頓

F2=0#另一端的力，單位：牛頓

#計算力向量

F=force_vector(F1,F2)

print(F)3.4桿單元的實例分析3.4.1問題描述考慮一根長度為1.0米的桿，截面積為0.01平方米，彈性模量為200吉帕。桿的一端固定，另一端受到1000牛頓的拉力。求桿的位移。3.4.2分析步驟計算剛度矩陣K。應(yīng)用力邊界條件，得到力向量F。解線性方程組Ku=F3.4.3代碼示例importnumpyasnp

#定義參數(shù)

A=0.01#截面積，單位：平方米

E=200e9#彈性模量，單位：帕斯卡

L=1.0#桿單元長度，單位：米

F1=1000#一端的力，單位：牛頓

F2=0#另一端的力，單位：牛頓

#計算剛度矩陣

K=stiffness_matrix(A,E,L)

#應(yīng)用力邊界條件

F=force_vector(F1,F2)

#解線性方程組

#由于一端固定，u1=0，因此我們只需要求解u2

#修改剛度矩陣和力向量

K_mod=[[K[1][1]]]

F_mod=[F[1]]

#使用numpy求解

u=np.linalg.solve(K_mod,F_mod)

print(u)在上述代碼中，我們首先計算了剛度矩陣K和力向量F。然后，考慮到一端固定，我們修改了剛度矩陣和力向量，只保留了與自由端相關(guān)的行和列。最后，使用numpy庫的linalg.solve函數(shù)求解線性方程組，得到自由端的位移u24梁單元分析4.1梁單元的數(shù)學模型梁單元在結(jié)構(gòu)力學中用于模擬承受橫向載荷的長條形結(jié)構(gòu)。其數(shù)學模型基于歐拉-伯努利梁理論，該理論假設(shè)梁是細長的，且橫向變形遠小于軸向變形。梁的變形由撓度函數(shù)描述，該函數(shù)表示梁在載荷作用下沿其長度方向的垂直位移。撓度函數(shù)滿足以下微分方程：E其中，E是彈性模量，I是截面慣性矩，w是撓度，qx4.2梁單元的剛度矩陣推導梁單元的剛度矩陣描述了梁單元在不同節(jié)點處的力與位移之間的關(guān)系。對于一個兩端固定的梁單元，其剛度矩陣為4x4矩陣，考慮兩個節(jié)點的轉(zhuǎn)角和撓度。剛度矩陣的推導基于虛功原理，即外力所做的虛功等于結(jié)構(gòu)內(nèi)部應(yīng)力所做的虛功。通過將梁單元的微分方程轉(zhuǎn)化為能量形式，可以得到梁單元的剛度矩陣。k其中，kij表示在節(jié)點i處施加單位力時，在節(jié)點4.3梁單元的轉(zhuǎn)角與撓度計算在有限元分析中，梁單元的轉(zhuǎn)角和撓度計算是通過求解結(jié)構(gòu)的整體剛度方程來實現(xiàn)的。整體剛度方程為：K其中，K是結(jié)構(gòu)的整體剛度矩陣，u是節(jié)點位移向量，F(xiàn)是節(jié)點力向量。通過將所有梁單元的剛度矩陣和載荷向量組合成整體矩陣和向量，可以求解出節(jié)點位移向量，進而得到梁單元的轉(zhuǎn)角和撓度。4.3.1示例代碼假設(shè)我們有一個簡單的梁單元，兩端固定，長度為L，彈性模量為E，截面慣性矩為I。我們使用Python和NumPy庫來計算該梁單元的剛度矩陣，并求解在端部施加的集中力作用下的轉(zhuǎn)角和撓度。importnumpyasnp

#定義梁單元的參數(shù)

E=200e9#彈性模量，單位：Pa

I=1e-4#截面慣性矩，單位：m^4

L=1.0#梁單元的長度，單位：m

#計算梁單元的剛度矩陣

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

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

[-12,-6*L,12,-6*L],

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

#定義節(jié)點力向量

F=np.array([0,-1000,0,0])#單位：N

#定義節(jié)點位移邊界條件

u=np.array([0,0,0,0])#單位：m

#求解節(jié)點位移向量

u=np.linalg.solve(k,F)

#輸出轉(zhuǎn)角和撓度

print("轉(zhuǎn)角：",u[0],u[2])

print("撓度：",u[1],u[3])4.3.2代碼解釋在上述代碼中，我們首先定義了梁單元的參數(shù)，包括彈性模量E、截面慣性矩I和長度L。然后，我們使用這些參數(shù)計算了梁單元的剛度矩陣。接下來，我們定義了節(jié)點力向量F，其中在第一個節(jié)點處施加了-1000N的集中力。由于梁的兩端固定，我們設(shè)定了節(jié)點位移邊界條件u為零。最后，我們使用NumPy的linalg.solve函數(shù)求解了節(jié)點位移向量，并輸出了轉(zhuǎn)角和撓度。4.4梁單元的實例分析4.4.1示例：懸臂梁分析考慮一個懸臂梁，長度為L，一端固定，另一端自由。在自由端施加一個集中力F。我們使用有限元法分析該梁的轉(zhuǎn)角和撓度。importnumpyasnp

#定義梁單元的參數(shù)

E=200e9#彈性模量，單位：Pa

I=1e-4#截面慣性矩，單位：m^4

L=1.0#梁單元的長度，單位：m

#計算梁單元的剛度矩陣

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

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

[-12,-6*L,12,-6*L],

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

#定義節(jié)點力向量

F=np.array([0,0,0,-1000])#單位：N

#定義節(jié)點位移邊界條件

u=np.array([0,0,0,0])#單位：m

#由于一端固定，另一端自由，我們只需要求解自由端的位移

#因此，我們修改剛度矩陣和力向量，只保留自由端的信息

k_free=k[2:,2:]

F_free=F[2:]

#求解節(jié)點位移向量

u_free=np.linalg.solve(k_free,F_free)

#輸出轉(zhuǎn)角和撓度

print("自由端轉(zhuǎn)角：",u_free[0])

print("自由端撓度：",u_free[1])4.4.2代碼解釋在這個懸臂梁的分析中，我們同樣定義了梁單元的參數(shù)，并計算了剛度矩陣。但是，由于懸臂梁的一端固定，另一端自由，我們只需要求解自由端的位移。因此，我們修改了剛度矩陣和力向量，只保留了自由端的信息。最后，我們求解了自由端的位移，并輸出了轉(zhuǎn)角和撓度。通過以上示例，我們可以看到，有限元法在梁單元分析中的應(yīng)用，包括剛度矩陣的計算和節(jié)點位移的求解。這些方法可以擴展到更復雜的結(jié)構(gòu)和載荷情況，為結(jié)構(gòu)力學的數(shù)值分析提供了強大的工具。5FEM在桿單元與梁單元中的應(yīng)用5.1桿單元與梁單元的組合結(jié)構(gòu)分析在結(jié)構(gòu)力學中，桿單元和梁單元是有限元分析(FEM)中最基本的單元類型。桿單元主要用于分析軸向受力的結(jié)構(gòu)，如拉桿、壓桿等，而梁單元則適用于分析受彎矩、剪力和扭矩的結(jié)構(gòu)，如橋梁、梁等。組合結(jié)構(gòu)分析是指將這兩種單元結(jié)合使用，以分析更復雜、更實際的結(jié)構(gòu)問題。5.1.1桿單元分析桿單元的分析基于一維桿件的力學原理，其主要考慮的是軸向力。在FEM中，桿單元的剛度矩陣是一個2x2的矩陣，表示兩端節(jié)點的位移與力之間的關(guān)系。例如，對于一個簡單的桿單元，其剛度矩陣可以表示為：#桿單元的剛度矩陣示例

#假設(shè)桿的彈性模量為E，截面積為A，長度為L

E=200e9#彈性模量，單位：Pa

A=0.01#截面積，單位：m^2

L=1.0#長度，單位：m

#桿單元的剛度矩陣

k=E*A/L*np.array([[1,-1],[-1,1]])5.1.2梁單元分析梁單元的分析則更為復雜，需要考慮彎矩、剪力和扭矩的影響。梁單元的剛度矩陣是一個6x6的矩陣，表示兩端節(jié)點的位移、轉(zhuǎn)角與力、彎矩之間的關(guān)系。例如，對于一個簡單的梁單元，其剛度矩陣可以表示為：#梁單元的剛度矩陣示例

#假設(shè)梁的彈性模量為E，慣性矩為I，長度為L

E=200e9#彈性模量，單位：Pa

I=0.001#慣性矩，單位：m^4

L=1.0#長度，單位：m

#梁單元的剛度矩陣

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

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

[-12,-6*L,12,-6*L],

[6*L,2*L**2,-6*L,4*L**2]])5.1.3組合結(jié)構(gòu)分析在實際應(yīng)用中，結(jié)構(gòu)往往由多種類型的單元組成，如桿單元和梁單元的組合。這種組合結(jié)構(gòu)的分析需要將所有單元的剛度矩陣組合成一個全局剛度矩陣，然后求解整個結(jié)構(gòu)的位移和內(nèi)力。例如，假設(shè)我們有一個由一個桿單元和一個梁單元組成的結(jié)構(gòu)，其全局剛度矩陣可以表示為：#全局剛度矩陣示例

#假設(shè)結(jié)構(gòu)由一個桿單元和一個梁單元組成

#桿單元的剛度矩陣為k1，梁單元的剛度矩陣為k2

k1=E*A/L*np.array([[1,-1],[-1,1]])

k2=E*I/L**3*np.array([[12,6*L,-12,6*L],

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

[-12,-6*L,12,-6*L],

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

#組合結(jié)構(gòu)的全局剛度矩陣

K=np.zeros((4,4))

K[:2,:2]=k1

K[2:,2:]=k25.2復雜結(jié)構(gòu)的FEM建模策略對于復雜結(jié)構(gòu)的FEM建模，需要考慮的因素更多，包括結(jié)構(gòu)的幾何形狀、材料屬性、邊界條件、載荷分布等。建模策略通常包括以下步驟：幾何建模：使用CAD軟件創(chuàng)建結(jié)構(gòu)的幾何模型。網(wǎng)格劃分：將幾何模型劃分為多個單元，單元的大小和形狀需要根據(jù)結(jié)構(gòu)的特性和分析的精度要求來確定。材料屬性和邊界條件：為每個單元指定材料屬性，如彈性模量、泊松比等，并定義邊界條件，如固定端、自由端等。載荷施加：在結(jié)構(gòu)上施加載荷，如重力、風力、地震力等。求解和后處理：求解結(jié)構(gòu)的位移和內(nèi)力，然后進行后處理，如繪制應(yīng)力云圖、變形圖等。5.3FEM軟件操作指南FEM軟件如ANSYS、ABAQUS、NASTRAN等，提供了強大的建模和分析工具，可以大大簡化復雜結(jié)構(gòu)的FEM分析過程。以下是一個使用Python的FEniCS庫進行FEM分析的基本操作指南：#導入必要的庫

fromfenicsimport*

#創(chuàng)建網(wǎng)格

mesh=UnitSquareMesh(8,8)

#定義函數(shù)空間

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

#定義邊界條件

defboundary(x,on_boundary):

returnon_boundary

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

#定義材料屬性和載荷

E=10.0

nu=0.3

f=Constant((0,-10))

#定義變分問題

u=TrialFunction(V)

v=TestFunction(V)

a=inner(nu*grad(u),grad(v))*dx

L=inner(f,v)*dx

#求解

u=Function(V)

solve(a==L,u,bc)這個例子中，我們創(chuàng)建了一個單位正方形的網(wǎng)格，定義了一個向量函數(shù)空間，然后定義了邊界條件、材料屬性和載荷。最后，我們定義了一個變分問題，并求解了結(jié)構(gòu)的位移。這只是一個非?；A(chǔ)的例子，實際的FEM分析會更復雜，需要根據(jù)具體的結(jié)構(gòu)和分析需求來調(diào)整。6非線性分析在桿單元與梁單元中的應(yīng)用非線性分析在結(jié)構(gòu)力學中是一個關(guān)鍵領(lǐng)域，尤其當結(jié)構(gòu)承受大變形、大應(yīng)變或材料行為超出線性范圍時。在桿單元與梁單元的分析中，非線性效應(yīng)可能源于幾何非線性、材料非線性或兩者結(jié)合。6.1幾何非線性幾何非線性分析考慮了結(jié)構(gòu)變形對自身幾何形狀的影響。在大變形情況下，結(jié)構(gòu)的初始幾何形狀與變形后的形狀差異顯著，不能忽略。例如，當梁單元承受較大的橫向載荷時，其軸向變形將影響橫向變形的計算，這種效應(yīng)稱為“P-δ”效應(yīng)。6.1.1示例：使用Python進行幾何非線性分析假設(shè)我們有一個簡單的梁單元，承受橫向載荷。我們將使用Python和一個流行的有限元分析庫FEniCS來模擬這一過程。fromfenicsimport*

#創(chuàng)建網(wǎng)格和函數(shù)空間

mesh=UnitIntervalMesh(10)

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

#定義邊界條件

defboundary(x,on_boundary):

returnon_boundary

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

#定義變量

u=TrialFunction(V)

v=TestFunction(V)

#定義材料屬性和載荷

E=1e3#彈性模量

nu=0.3#泊松比

rho=1.0#密度

f=Constant((0,-10))#橫向載荷

#定義應(yīng)變和應(yīng)力

defepsilon(u):

returnsym(nabla_grad(u))

defsigma(u):

returnE/(1+nu)*epsilon(u)

#定義弱形式

F=inner(sigma(u),epsilon(v))*dx-inner(f,v)*dx

#解非線性問題

solve(F==0,u,bc)

#可視化結(jié)果

plot(u)

interactive()在這個例子中，我們首先創(chuàng)建了一個單元區(qū)間網(wǎng)格，并定義了向量函數(shù)空間。然后，我們設(shè)置了邊界條件，確保梁的兩端固定。通過定義材料屬性和載荷，我們能夠計算梁在非線性條件下的變形。最后，我們使用FEniCS的solve函數(shù)來求解非線性方程，并可視化結(jié)果。6.2材料非線性材料非線性分析考慮了材料在不同應(yīng)力水平下的行為變化。例如，當應(yīng)力超過材料的屈服點時，材料將進入塑性狀態(tài)，其應(yīng)力-應(yīng)變關(guān)系不再是線性的。在桿單元和梁單元的分析中，這種效應(yīng)尤其重要，因為它影響結(jié)構(gòu)的承載能力和安全性。6.2.1示例：使用Python進行材料非線性分析我們將使用FEniCS來模擬一個梁單元在材料進入塑性狀態(tài)時的行為。fromfenicsimport*

#創(chuàng)建網(wǎng)格和函數(shù)空間

mesh=UnitIntervalMesh(10)

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

#定義邊界條件

defboundary(x,on_boundary):

returnon_boundary

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

#定義變量

u=TrialFunction(V)

v=TestFunction(V)

#定義材料屬性和載荷

E=1e3#彈性模量

nu=0.3#泊松比

sigma_y=100#屈服應(yīng)力

rho=1.0#密度

f=Constant((0,-10))#橫向載荷

#定義應(yīng)變和應(yīng)力

defepsilon(u):

returnsym(nabla_grad(u))

defsigma(u):

e=epsilon(u)

returnE/(1+nu)*e-E/(1+nu)*sigma_y*(e>sigma_y/E)

#定義弱形式

F=inner(sigma(u),epsilon(v))*dx-inner(f,v)*dx

#解非線性問題

solve(F==0,u,bc)

#可視化結(jié)果

plot(u)

interactive()在這個例子中，我們定義了一個塑性應(yīng)力-應(yīng)變關(guān)系，當應(yīng)變超過屈服應(yīng)變時，應(yīng)力將不再線性增加。通過FEniCS的solve函數(shù)，我們能夠求解梁單元在材料非線性條件下的變形。7材料屬性對結(jié)構(gòu)響應(yīng)的影響材料屬性，如彈性模量、泊松比和屈服強度，對結(jié)構(gòu)的響應(yīng)有著直接的影響。在有限元分析中，這些屬性的準確設(shè)定對于預測結(jié)構(gòu)行為至關(guān)重要。7.1彈性模量彈性模量（E）是材料在彈性范圍內(nèi)抵抗變形的能力的度量。較高的彈性模量意味著材料更硬，變形更小。7.1.1示例：彈性模量對梁單元變形的影響我們使用FEniCS來模擬不同彈性模量下的梁單元變形。fromfenicsimport*

#創(chuàng)建網(wǎng)格和函數(shù)空間

mesh=UnitIntervalMesh(10)

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

#定義邊界條件

defboundary(x,on_boundary):

returnon_boundary

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

#定義變量

u=TrialFunction(V)

v=TestFunction(V)

#定義材料屬性和載荷

nu=0.3#泊松比

rho=1.0#密度

f=Constant((0,-10))#橫向載荷

#定義應(yīng)變和應(yīng)力

defepsilon(u):

returnsym(nabla_grad(u))

defsigma(u,E):

returnE/(1+nu)*epsilon(u)

#定義弱形式

F=inner(sigma(u,1e3),epsilon(v))*dx-inner(f,v)*dx

#解非線性問題

u1=Function(V)

solve(F==0,u1,bc)

#改變彈性模量

F=inner(sigma(u,1e4),epsilon(v))*dx-inner(f,v)*dx

#解非線性問題

u2=Function(V)

solve(F==0,u2,bc)

#可視化結(jié)果

plot(u1,title='E=1e3')

plot(u2,title='E=1e4')

interactive()在這個例子中，我們首先求解了彈性模量為1e3的梁單元變形，然后將彈性模量提高到1e4，再次求解。通過比較兩個結(jié)果，我們可以直觀地看到彈性模量對梁單元變形的影響。7.2泊松比泊松比（ν）描述了材料在拉伸或壓縮時橫向變形與縱向變形的比例。不同的泊松比會影響結(jié)構(gòu)的橫向變形。7.2.1示例：泊松比對梁單元橫向變形的影響我們使用FEniCS來模擬不同泊松比下的梁單元橫向變形。fromfenicsimport*

#創(chuàng)建網(wǎng)格和函數(shù)空間

mesh=UnitIntervalMesh(10)

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

#定義邊界條件

defboundary(x,on_boundary):

returnon_boundary

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

#定義變量

u=TrialFunction(V)

v=TestFunction(V)

#定義材料屬性和載荷

E=1e3#彈性模量

rho=1.0#密度

f=Constant((0,-10))#橫向載荷

#定義應(yīng)變和應(yīng)力

defepsilon(u):

returnsym(nabla_grad(u))

defsigma(u,nu):

returnE/(1+nu)*epsilon(u)

#定義弱形式

F=inner(sigma(u,0.3),epsilon(v))*dx-inner(f,v)*dx

#解非線性問題

u1=Function(V)

solve(F==0,u1,bc)

#改變泊松比

F=inner(sigma(u,0.4),epsilon(v))*dx-inner(f,v)*dx

#解非線性問題

u2=Function(V)

solve(F==0,u2,bc)

#可視化結(jié)果

plot(u1,title='ν=0.3')

plot(u2,title='ν=0.4')

interactive()在這個例子中，我們首先求解了泊松比為0.3的梁單元變形，然后將泊松比提高到0.4，再次求解。通過比較兩個結(jié)果，我們可以看到泊松比對梁單元橫向變形的影響。8FEM結(jié)果的后處理與解釋有限元分析的結(jié)果需要通過后處理來解釋和可視化，以便于理解和分析。這包括應(yīng)力、應(yīng)變、位移等結(jié)果的可視化，以及對結(jié)果的進一步分析，如安全系數(shù)計算。8.1示例：使用Python進行FEM結(jié)果的后處理我們將使用FEniCS來模擬一個梁單元的變形，并使用matplotlib庫來可視化結(jié)果。importmatplotlib.pyplotasplt

fromfenicsimport*

#創(chuàng)建網(wǎng)格和函數(shù)空間

mesh=UnitIntervalMesh(10)

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

#定義邊界條件

defboundary(x,on_boundary):

returnon_boundary

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

#定義變量

u=TrialFunction(V)

v=TestFunction(V)

#定義材料屬性和載荷

E=1e3#彈性模量

nu=0.3#泊松比

rho=1.0#密度

f=Constant((0,-10))#橫向載荷

#定義應(yīng)變和應(yīng)力

defepsilon(u):

returnsym(nabla_grad(u))

defsigma(u):

returnE/(1+nu)*epsilon(u)

#定義弱形式

F=inner(sigma(u),epsilon(v))*dx-inner(f,v)*dx

#解非線性問題

u_solution=Function(V)

solve(F==0,u_solution,bc)

#可視化結(jié)果

plt.figure()

plot(u_solution)

plt.title('梁單元變形')

plt.xlabel('長度')

plt.ylabel('位移')

plt.show()在這個例子中，我們使用matplotlib庫來可視化梁單元的變形結(jié)果。通過繪制位移圖，我們可以直觀地看到梁單元在載荷作用下的變形情況。此外，我們還可以進一步分析應(yīng)力、應(yīng)變等結(jié)果，以評估結(jié)構(gòu)的安全性和性能。9案例研究9.1橋梁結(jié)構(gòu)的FEM分析9.1.1原理與內(nèi)容橋梁結(jié)構(gòu)的有限元分析(FEM)是一種數(shù)值方法，用于預測橋梁在各種載荷條件下的行為。它將橋梁結(jié)構(gòu)分解為多個小的、簡單的部分，即單元，每個單元的行為可以通過已知的物理定律來描述。這些單元通過節(jié)點連接，形成整個結(jié)構(gòu)的模型。通過求解單元的平衡方程，可以得到整個結(jié)構(gòu)的響應(yīng)，包括位移、應(yīng)力和應(yīng)變。示例：橋梁的靜態(tài)分析假設(shè)我們有一個簡單的橋梁模型，由多個桿單元和梁單元組成。我們將使用Python的numpy庫和scipy庫來執(zhí)行靜態(tài)分析。importnumpyasnp

fromscipy.sparseimportlil_matrix

fromscipy.sparse.linalgimportspsolve

#定義結(jié)構(gòu)參數(shù)

n_nodes=10#節(jié)點數(shù)量

n_elements=9#元素數(shù)量

E=200e9#彈性模量，單位：Pa

A=0.1#截面積，單位：m^2

I=0.001#慣性矩，單位：m^4

L=10#橋梁總長度，單位：m

node_spacing=L/(n_nodes-1)#節(jié)點間距

#創(chuàng)建節(jié)點坐標

nodes=np.linspace(0,L,n_nodes)

#創(chuàng)建單元連接矩陣

elements=np.zeros((n_elements,2),dtype=int)

foriinrange(n_elements):

elements[i]=[i,i+1]

#創(chuàng)建剛度矩陣

K=lil_matrix((2*n_nodes,2*n_nodes))

fori,(n1,n2)inenumerate(elements):

#桿單元剛度矩陣

k_rod=(E*A/node_spacing)*np.array([[1,-1],[-1,1]])

#梁單元剛度矩陣

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

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

[-12,-6*node_spacing,12,-6*node_spacing],

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

#將局部剛度矩陣轉(zhuǎn)換為全局坐標系

#注意：這里簡化了轉(zhuǎn)換過程，實際應(yīng)用中需要考慮單元的方向和坐標系轉(zhuǎn)換

K[2*n1:2*n1+2,2*n1:2*n1+2]+=k_rod

K[2*n1:2*n1+4,2*n1:2*n1+4]+=k_beam

#應(yīng)用邊界條件

#假設(shè)第一個節(jié)點固定，最后一個節(jié)點有向下的力

K[0,:]=0

K[0,0]=1

K[-1,:]=0

K[-1,-1]=1

#定義載荷向量

F=np.zeros(2*n_nodes)

F[-1]=-1000#單位：N

#求解位移向量

U=spsolve(K.tocsr(),F)

#計算應(yīng)力和應(yīng)變

#這里簡化了計算過程，實際應(yīng)用中需要根據(jù)位移向量和單元類型計算應(yīng)力和應(yīng)變

stresses=np.zeros(n_elements)

fori,(n1,n2)inenumerate(elements):

#桿單元應(yīng)力計算

stresses[i]=E*(U[2*n2]-U[2*n1])/node_spacing

#輸出結(jié)果

print("節(jié)點位移：",U)

print("單元應(yīng)力：",stresses)9.1.2描述上述代碼示例展示了如何使用有限元方法分析一個由桿單元和梁單元組成的橋梁結(jié)構(gòu)。首先，我們定義了結(jié)構(gòu)的基本參數(shù)，如節(jié)點數(shù)量、元素數(shù)量、彈性模量、截面積、慣性矩和橋梁總長度。然后，我們創(chuàng)建了節(jié)點坐標和元素連接矩陣，用于描述結(jié)構(gòu)的幾何形狀。接下來，我們創(chuàng)建了剛度矩陣K，它是一個稀疏矩陣，用于存儲結(jié)構(gòu)的剛度信息。對于每個元素，我們計算了桿單元和梁單元的局部剛度矩陣，并將它們轉(zhuǎn)換到全局坐標系中，然后累加到K中。這里我們簡化了坐標系轉(zhuǎn)換的過程，實際應(yīng)用中需要考慮單元的方向和坐標系轉(zhuǎn)換。我們應(yīng)用了邊界條件，假設(shè)第一個節(jié)點固定，最后一個節(jié)點有向下的力。然后，我們定義了載荷向量F，并使用scipy.sparse.linalg.spsolve函數(shù)求解位移向量U。最后，我們計算了每個單元的應(yīng)力，這里也簡化了計算過程，實際應(yīng)用中需要根據(jù)位移向量和單元類型計算應(yīng)力和應(yīng)變。我們輸出了節(jié)點位移和單元應(yīng)力的結(jié)果。9.2高層建筑框架的桿單元與梁單元分析9.2.1原理與內(nèi)容高層建筑框架的有限元分析通常涉及復雜的三維結(jié)構(gòu)，包括多個桿單元和梁單元。分析的目標是評估結(jié)構(gòu)在地震、風載荷等動態(tài)載荷下的響應(yīng)，以及在重力載荷下的靜態(tài)響應(yīng)。這種分析需要考慮結(jié)構(gòu)的幾何非線性、材料非線性和接觸非線性。示例：高層建筑框架的靜態(tài)分析假設(shè)我們有一個簡單的高層建筑框架模型，由多個桿單元和梁單元組成。我們將使用Python的numpy庫和scipy庫來執(zhí)行靜態(tài)分析。importnumpyasnp

fromscipy.sparseimportlil_matrix

fromscipy.sparse.linalgimportspsolve

#定義結(jié)構(gòu)參數(shù)

n_nodes=20#節(jié)點數(shù)量

n_elements=18#元素數(shù)量

E=30e9#彈性模量，單位：Pa

A=0.2#截面積，單位：m^2

I=0.002#慣性矩，單位：m^4

H=100#建筑總高度，單位：m

node_spacing=H/(n_nodes-1)#節(jié)點間距

#創(chuàng)建節(jié)點坐標

nodes=np.linspace(0,H,n_nodes)

#創(chuàng)建單元連接矩陣

elements=np.zeros((n_elements,2),dtype=int)

foriinrange(n_elements):

elements[i]=[i,i+1]

#創(chuàng)建剛度矩陣

K=lil_matrix((3*n_nodes,3*n_nodes))

fori,(n1,n2)inenumerate(elements):

#桿單元剛度矩陣

k_rod=(E*A/node_spacing)*np.array([[1,-1],[-1,1]])

#梁單元剛度矩陣

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

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

[-12,-6*node_spacing,12,-6*node_spacing],

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

#將局部剛度矩陣轉(zhuǎn)換為全局坐標系

#注意：這里簡化了轉(zhuǎn)換過程，實際應(yīng)用中需要考慮單元的方向和坐標系轉(zhuǎn)換

K[3*n1:3*n1+2,3*n1:3*n1+2]+=k_rod

K[3*n1:3*n1+4,3*n1:3*n1+4]+=k_beam

#應(yīng)用邊界條件

#假設(shè)第一個節(jié)點固定，最后一個節(jié)點有向下的力

K[0,:]=0

K[0,0]=1

K[1,:]=0

K[1,1]=1

K[-1,:]=0

K[-1,-1]=1

#定義載荷向量

F=np.zeros(3*n_nodes)

F[-1]=-10000#單位：N

#求解位移向量

U=spsolve(K.tocsr(),F)

#計算應(yīng)力和應(yīng)變

#這里簡化了計算過程，實際應(yīng)用中需要根據(jù)位移向量和單元類型計算應(yīng)力和應(yīng)變

stresses=np.zeros(n_elements)

fori,(n1,n2)inenumerate(elements):

#桿單元應(yīng)力計算

stresses[i]=E*(U[3*n2]-U[3*n1])/node_spacing

#輸出結(jié)果

print("節(jié)點位移：",U)

print("單元應(yīng)力：",stresses)9.2.2描述上述代碼示例展示了如何使用有限元方法分析一個由桿單元和梁單元組成的高層建筑框架結(jié)構(gòu)。與橋梁結(jié)構(gòu)分析類似，我們首先定義了結(jié)構(gòu)的基本參數(shù)，如節(jié)點數(shù)量、元素數(shù)量、彈性模量、截面積、慣性矩和建筑總高度。然后，我們創(chuàng)建了節(jié)點坐標和元素連接矩陣，用于描述結(jié)構(gòu)的幾何形狀。接下來，我們創(chuàng)建了剛度矩陣K，它是一個稀疏矩陣，用于存儲結(jié)構(gòu)的剛度信息。對于每個元素，我們計算了桿單元和梁單元的局部剛度矩陣，并將它們轉(zhuǎn)換到全局坐標系中，然后累加到K中。這里我們簡化了坐標系轉(zhuǎn)換的過程，實際應(yīng)用中需要考慮單元的方向和坐標系轉(zhuǎn)換。我們應(yīng)用了邊界條件，假設(shè)第一個節(jié)點在X和Y方向固定，最后一個節(jié)點有向下的力。然后，我們定義了載荷向量F，并使用scipy.sparse.linalg.spsolve函數(shù)求解位移向量U。最后，我們計算了每個單元的應(yīng)力，這里也簡化了計算過程，實際應(yīng)用中需要根據(jù)位移向量和單元類型計算應(yīng)力和應(yīng)變。我們輸出了節(jié)點位移和單元應(yīng)力的結(jié)果。9.3機械臂的動態(tài)響應(yīng)分析9.3.1原理與內(nèi)容機械臂的動態(tài)響應(yīng)分析通常涉及復雜的運動學和動力學問題。分析的目標是評估機械臂在運動過程中的響應(yīng)，包括位移、速度、加速度和應(yīng)力。這種分析需要考慮結(jié)構(gòu)的慣性、阻尼和外部載荷的影響。示例：機械臂的動態(tài)響應(yīng)分析假設(shè)我們有一個簡單的機械臂模型，由多個桿單元和梁單元組成。我們將使用Python的numpy庫和scipy庫來執(zhí)行動態(tài)響應(yīng)分析。importnumpyasnp

fromscipy.sparseimportlil_matrix

fromscipy.sparse.linalgimportspsolve

fromegrateimportsolve_ivp

#定義結(jié)構(gòu)參數(shù)

n_nodes=10#節(jié)點數(shù)量

n_elements=9#元素數(shù)量

E=200e9#彈性模量，單位：Pa

A=0.1#截面積，單位：m^2

I=0.001#慣性矩，單位：m^4

L=10#機械臂總長度，單位：m

node_spacing=L/(n_nodes-1)#節(jié)點間距

rho=7850#材料密度，單位：kg/m^3

m=rho*A*node_spacing#單元質(zhì)量，單位：kg

#創(chuàng)建節(jié)點坐標

nodes=np.linspace(0,L,n_nodes)

#創(chuàng)建單元連接矩陣

elements=np.zeros((n_elements,2),dtype=int)

foriinrange(n_elements):

elements[i]=[i,i+1]

#創(chuàng)建剛度矩陣和質(zhì)量矩陣

K=lil_matrix((2*n_nodes,2*n_nodes))

M=lil_matrix((2*n_nodes,2*n_nodes))

fori,(n1,n2)inenumerate(elements):

#桿單元剛度矩陣

k_rod=(E*A/node_spacing)*np.array([[1,-1],[-1,1]])

#梁單元剛度矩陣

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

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

[-12,-6*node_spacing,12,-6*node_spacing],

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

#桿單元質(zhì)量矩陣

m_rod=m*np.array([[1,0],[0,1]])

#梁單元質(zhì)量矩陣

m_beam=m*np.array([[156,22*node_spacing,54,-13*node_spacing],

[22*node_spacing,4*node_spacing**2,-22*node_spacing,6*node_spacing**2],

[54,-22*node_spacing,156,-22*node_spacing],

[-13*node_spacing,6*node_spacing**2,-22*node_spacing,4*node_spacing**2]])

#將局部剛度矩陣和質(zhì)量矩陣轉(zhuǎn)換為全局坐標系

#注意：這里簡化了轉(zhuǎn)換過程，實際應(yīng)用中需要考慮單元的方向和坐標系轉(zhuǎn)換

K[2*n1:2*n1+2,2*n1:2*n1+2]+=k_rod

K[2*n1:2*n1+4,2*n1:2*n1+4]+=k_beam

M[2*n1:2*n1+2,2*n1:2*n1+2]+=m_rod

M[2*n1:2*n1+4,2*n1:2*n1+4]+=m_beam

#應(yīng)用邊界條件

#假設(shè)第一個節(jié)點固定

K[0,:]=0

K[0,0]=1

M[0,:]=0

M[0,0]=1

#定義動態(tài)方程

defdynamic_eq(t,y):

u=y[:2*n_nodes]

v=y[2*n_nodes:]

du=np.zeros(2*n_nodes)

dv=np.zeros(2*n_nodes)

#應(yīng)用外部載荷

F=np.zeros(2*n_nodes)

F[-1]=-1000*np.sin(2*np.pi*t)#單位：N

#求解加速度

dv=spsolve(M.tocsr(),F-K.dot(u)-M.dot(v))

#更新速度和位移

du=v

returnnp.concatenate((du,