彈性力學(xué)數(shù)值方法：有限元法(FEM)：彈性力學(xué)基礎(chǔ)理論

彈性力學(xué)數(shù)值方法：有限元法(FEM)：彈性力學(xué)基礎(chǔ)理論1彈性力學(xué)基礎(chǔ)1.1應(yīng)力與應(yīng)變1.1.1應(yīng)力張量的概念在彈性力學(xué)中，應(yīng)力張量描述了物體內(nèi)部各點(diǎn)處的應(yīng)力狀態(tài)。它是一個(gè)二階張量，可以表示為一個(gè)3x3的矩陣，其中包含了正應(yīng)力和剪應(yīng)力的分量。正應(yīng)力表示垂直于截面的力，而剪應(yīng)力表示平行于截面的力。應(yīng)力張量的分量可以通過(guò)以下公式計(jì)算：σ其中，σij是應(yīng)力張量的分量，F(xiàn)i是在i方向上的力，A1.1.2應(yīng)變張量的概念應(yīng)變張量描述了物體在受力作用下的變形程度。它同樣是一個(gè)二階張量，可以表示為3x3的矩陣，包含了線應(yīng)變和剪應(yīng)變的分量。線應(yīng)變表示物體在某一方向上的長(zhǎng)度變化，而剪應(yīng)變表示物體形狀的改變。應(yīng)變張量的分量可以通過(guò)以下公式計(jì)算：?其中，?ij是應(yīng)變張量的分量，ui和uj是位移分量，1.1.3胡克定律與材料屬性胡克定律描述了線性彈性材料的應(yīng)力與應(yīng)變之間的關(guān)系。對(duì)于各向同性材料，胡克定律可以表示為：σ其中，Ciσ其中，λ和μ是拉梅常數(shù)，δi1.2平衡方程與邊界條件1.2.1彈性體的平衡方程平衡方程描述了物體內(nèi)部的力平衡條件。在彈性力學(xué)中，平衡方程可以表示為：?其中，σij是應(yīng)力張量的分量，fi是體積力的分量，ρ1.2.2位移邊界條件與應(yīng)力邊界條件在有限元分析中，邊界條件是解決問(wèn)題的關(guān)鍵。位移邊界條件指定物體在邊界上的位移，而應(yīng)力邊界條件指定物體在邊界上的應(yīng)力。例如，如果物體的一端被固定，那么在該端的位移邊界條件可以表示為：u如果物體的一端受到外力作用，那么在該端的應(yīng)力邊界條件可以表示為：σ其中，nj是邊界法向量的分量，t1.3能量原理1.3.1應(yīng)變能與外力功應(yīng)變能是物體在受力作用下變形時(shí)儲(chǔ)存的能量，可以表示為：U其中，V是物體的體積。外力功是外力對(duì)物體所做的功，可以表示為：W其中，?V是物體的邊界，d1.3.2最小勢(shì)能原理最小勢(shì)能原理指出，當(dāng)物體處于平衡狀態(tài)時(shí)，其勢(shì)能（應(yīng)變能與外力功之和）達(dá)到最小值。這一原理在有限元分析中被廣泛使用，用于求解物體的位移。1.3.3虛功原理虛功原理是彈性力學(xué)中的一個(gè)重要原理，它指出，對(duì)于任何虛位移，虛功的總和為零。這一原理可以用于驗(yàn)證有限元分析的正確性。1.4示例：計(jì)算彈性體的應(yīng)變能假設(shè)我們有一個(gè)長(zhǎng)方體彈性體，其尺寸為1m×1m×1mimportnumpyasnp

#定義材料屬性

E=1e5#彈性模量

nu=0.3#泊松比

mu=E/(2*(1+nu))#切變模量

lmbda=E*nu/((1+nu)*(1-2*nu))#拉梅常數(shù)

#定義應(yīng)力張量

sigma=np.array([[1e3,0,0],[0,0,0],[0,0,0]])

#計(jì)算應(yīng)變張量

epsilon=np.linalg.inv(2*mu*np.eye(3)+lmbda*np.ones((3,3)))@sigma

#定義體積

V=1*1*1

#計(jì)算應(yīng)變能

U=0.5*V*np.trace(sigma@epsilon)

print("應(yīng)變能為：",U)在這個(gè)例子中，我們首先定義了材料的彈性模量和泊松比，然后計(jì)算了切變模量和拉梅常數(shù)。接著，我們定義了應(yīng)力張量，并使用胡克定律計(jì)算了應(yīng)變張量。最后，我們定義了體積，并使用應(yīng)變能的公式計(jì)算了應(yīng)變能。1.5結(jié)論通過(guò)以上內(nèi)容，我們了解了彈性力學(xué)的基礎(chǔ)理論，包括應(yīng)力與應(yīng)變的概念、平衡方程與邊界條件、能量原理等。這些理論是有限元分析的基礎(chǔ)，對(duì)于理解和解決彈性力學(xué)問(wèn)題至關(guān)重要。2有限元法原理2.1離散化過(guò)程2.1.1網(wǎng)格劃分與節(jié)點(diǎn)編號(hào)在有限元分析中，首先需要將連續(xù)的結(jié)構(gòu)體離散化為一系列的單元。這個(gè)過(guò)程稱為網(wǎng)格劃分。每個(gè)單元由節(jié)點(diǎn)連接，而節(jié)點(diǎn)是計(jì)算位移、應(yīng)變和應(yīng)力的基本點(diǎn)。節(jié)點(diǎn)編號(hào)是為每個(gè)節(jié)點(diǎn)分配一個(gè)唯一的數(shù)字標(biāo)識(shí)，以便在計(jì)算過(guò)程中引用。例如，考慮一個(gè)簡(jiǎn)單的梁結(jié)構(gòu)，可以將其離散化為多個(gè)線性單元，每個(gè)單元有兩個(gè)節(jié)點(diǎn)。節(jié)點(diǎn)編號(hào)可以按順序從左至右進(jìn)行，確保每個(gè)節(jié)點(diǎn)都有一個(gè)唯一的編號(hào)。2.1.2單元與節(jié)點(diǎn)位移每個(gè)單元的位移由其節(jié)點(diǎn)位移決定。在二維問(wèn)題中，每個(gè)節(jié)點(diǎn)有兩個(gè)位移分量（x和y方向）。這些位移分量是有限元分析中的未知變量，通過(guò)求解線性方程組來(lái)確定。2.1.3應(yīng)變與應(yīng)力的插值應(yīng)變和應(yīng)力在單元內(nèi)部通過(guò)位移的插值函數(shù)來(lái)計(jì)算。對(duì)于線性單元，位移插值函數(shù)通常是線性的，這意味著應(yīng)變和應(yīng)力在單元內(nèi)部也是線性的。插值函數(shù)將節(jié)點(diǎn)位移映射到單元內(nèi)部的位移，從而計(jì)算出應(yīng)變和應(yīng)力。2.2剛度矩陣與載荷向量2.2.1局部剛度矩陣的推導(dǎo)局部剛度矩陣描述了單個(gè)單元內(nèi)部的力與位移之間的關(guān)系。它通過(guò)將單元的應(yīng)變能表示為位移的函數(shù)來(lái)推導(dǎo)。對(duì)于一個(gè)線性彈性單元，局部剛度矩陣可以通過(guò)以下公式計(jì)算：K其中，Ke是局部剛度矩陣，B是應(yīng)變位移矩陣，D是彈性矩陣，V2.2.2整體剛度矩陣的組裝整體剛度矩陣是通過(guò)將所有局部剛度矩陣組裝在一起得到的。這個(gè)過(guò)程涉及到將局部坐標(biāo)系下的剛度矩陣轉(zhuǎn)換到全局坐標(biāo)系，并將它們加到相應(yīng)的位置上。整體剛度矩陣描述了整個(gè)結(jié)構(gòu)的力與位移之間的關(guān)系。2.2.3等效節(jié)點(diǎn)載荷向量等效節(jié)點(diǎn)載荷向量是將作用在結(jié)構(gòu)上的外力轉(zhuǎn)換為作用在節(jié)點(diǎn)上的力。這個(gè)過(guò)程涉及到將外力通過(guò)單元的應(yīng)變位移矩陣和彈性矩陣轉(zhuǎn)換到節(jié)點(diǎn)載荷向量上。2.3求解過(guò)程2.3.1線性方程組的建立有限元分析最終會(huì)得到一個(gè)線性方程組，形式為：K其中，K是整體剛度矩陣，U是節(jié)點(diǎn)位移向量，F(xiàn)是等效節(jié)點(diǎn)載荷向量。2.3.2求解位移向量線性方程組可以通過(guò)多種數(shù)值方法求解，如直接求解法（如高斯消元法）或迭代求解法（如共軛梯度法）。求解得到的位移向量U包含了所有節(jié)點(diǎn)的位移信息。2.3.3計(jì)算應(yīng)變與應(yīng)力一旦得到了位移向量U，就可以通過(guò)應(yīng)變位移矩陣B和彈性矩陣D來(lái)計(jì)算每個(gè)單元的應(yīng)變和應(yīng)力：εσ其中，ε是應(yīng)變向量，σ是應(yīng)力向量。2.4示例：二維線性單元的有限元分析假設(shè)我們有一個(gè)簡(jiǎn)單的二維梁結(jié)構(gòu)，由兩個(gè)線性單元組成，每個(gè)單元有兩個(gè)節(jié)點(diǎn)。我們將使用Python和NumPy庫(kù)來(lái)演示如何進(jìn)行有限元分析。importnumpyasnp

#定義單元的節(jié)點(diǎn)編號(hào)

nodes_per_element=2

elements=np.array([[1,2],[2,3]])

#定義節(jié)點(diǎn)坐標(biāo)

node_coords=np.array([[0,0],[1,0],[2,0]])

#定義彈性矩陣（對(duì)于線性彈性材料）

D=np.array([[1,0],[0,1]])*100#假設(shè)彈性模量為100

#定義等效節(jié)點(diǎn)載荷向量

F=np.array([0,0,0,100,0,0])#假設(shè)在節(jié)點(diǎn)2的y方向有100的力

#定義局部剛度矩陣的推導(dǎo)函數(shù)

deflocal_stiffness_matrix(node1,node2,D):

L=np.sqrt((node2[0]-node1[0])**2+(node2[1]-node1[1])**2)

B=np.array([[-1/L,0,1/L,0],

[0,-1/L,0,1/L]])

Ke=np.dot(np.dot(B.T,D),B)*L

returnKe

#計(jì)算整體剛度矩陣

K=np.zeros((nodes_per_element*len(node_coords),nodes_per_element*len(node_coords)))

fore,elementinenumerate(elements):

Ke=local_stiffness_matrix(node_coords[element[0]-1],node_coords[element[1]-1],D)

foriinrange(nodes_per_element):

forjinrange(nodes_per_element):

K[(element[i]-1)*nodes_per_element:(element[i]-1)*nodes_per_element+nodes_per_element,

(element[j]-1)*nodes_per_element:(element[j]-1)*nodes_per_element+nodes_per_element]+=Ke[i*nodes_per_element:(i+1)*nodes_per_element,

j*nodes_per_element:(j+1)*nodes_per_element]

#求解位移向量

#假設(shè)節(jié)點(diǎn)1和節(jié)點(diǎn)3在x方向固定，因此位移為0

U=np.linalg.solve(K,F)

U[0]=0#節(jié)點(diǎn)1的x位移

U[2]=0#節(jié)點(diǎn)1的y位移

U[4]=0#節(jié)點(diǎn)3的x位移

U[5]=0#節(jié)點(diǎn)3的y位移

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

fore,elementinenumerate(elements):

B=np.array([[-1,0,1,0],

[0,-1,0,1]])

ue=U[(element[0]-1)*nodes_per_element:(element[1]-1)*nodes_per_element+nodes_per_element]

epsilon=np.dot(B,ue)

sigma=np.dot(D,epsilon)

print(f"單元{e+1}的應(yīng)變：{epsilon}")

print(f"單元{e+1}的應(yīng)力：{sigma}")在這個(gè)例子中，我們首先定義了單元的節(jié)點(diǎn)編號(hào)、節(jié)點(diǎn)坐標(biāo)、彈性矩陣和等效節(jié)點(diǎn)載荷向量。然后，我們定義了一個(gè)函數(shù)來(lái)計(jì)算局部剛度矩陣，并使用這個(gè)函數(shù)來(lái)組裝整體剛度矩陣。最后，我們求解了線性方程組得到位移向量，并計(jì)算了每個(gè)單元的應(yīng)變和應(yīng)力。請(qǐng)注意，這個(gè)例子是一個(gè)簡(jiǎn)化的二維線性彈性問(wèn)題，實(shí)際的有限元分析可能涉及到更復(fù)雜的單元類型、材料屬性和邊界條件。3有限元法應(yīng)用3.1平面應(yīng)力與平面應(yīng)變3.1.1平面應(yīng)力問(wèn)題的有限元分析平面應(yīng)力問(wèn)題通常出現(xiàn)在薄板結(jié)構(gòu)中，其中厚度方向的應(yīng)力可以忽略。在有限元分析中，我們使用平面應(yīng)力單元來(lái)模擬這類問(wèn)題。平面應(yīng)力單元假設(shè)應(yīng)力在厚度方向上為零，因此只有平面內(nèi)的應(yīng)力和應(yīng)變需要計(jì)算。3.1.1.1示例：平面應(yīng)力問(wèn)題的有限元分析假設(shè)我們有一個(gè)矩形薄板，其尺寸為100mmx50mm，厚度為1mm，材料為鋼，彈性模量為200GPa，泊松比為0.3。薄板的一端固定，另一端受到100N的水平拉力。我們將使用平面應(yīng)力單元來(lái)分析薄板的變形和應(yīng)力分布。#導(dǎo)入必要的庫(kù)

importnumpyasnp

fromfenicsimport*

#定義幾何參數(shù)

length=100.0

height=50.0

thickness=1.0

E=200e9#彈性模量

nu=0.3#泊松比

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

mesh=RectangleMesh(Point(0,0),Point(length,height),10,5)

#定義邊界條件

defleft_boundary(x,on_boundary):

returnnear(x[0],0.0)

defright_boundary(x,on_boundary):

returnnear(x[0],length)

#定義材料屬性

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

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

#定義變分形式

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

du=TrialFunction(V)

u_=TestFunction(V)

u=Function(V)

#應(yīng)力張量

defsigma(v):

returnlmbda*tr(eps(v))*Identity(2)+2*mu*eps(v)

#應(yīng)變張量

defeps(v):

returnsym(nabla_grad(v))

#定義邊界條件

bc_left=DirichletBC(V,Constant((0,0)),left_boundary)

bc_right=DirichletBC(V.sub(0),Constant(100/length),right_boundary)

#定義外力

f=Constant((0,0))

T=Constant((100,0))

#定義變分問(wèn)題

F=inner(sigma(du),eps(u_))*dx-inner(f,u_)*dx-inner(T,u_)*ds

#求解

solve(F==0,u,[bc_left,bc_right])

#輸出結(jié)果

file=File("displacement.pvd")

file<<u3.1.2平面應(yīng)變問(wèn)題的有限元分析平面應(yīng)變問(wèn)題適用于厚度遠(yuǎn)大于平面尺寸的結(jié)構(gòu)，如水壩或隧道壁。在這種情況下，厚度方向的應(yīng)變可以認(rèn)為是零，但應(yīng)力在厚度方向上可能存在。平面應(yīng)變單元考慮了這一特性，使得分析更加準(zhǔn)確。3.1.2.1示例：平面應(yīng)變問(wèn)題的有限元分析考慮一個(gè)水壩，其高度為50m，寬度為10m，材料為混凝土，彈性模量為30GPa，泊松比為0.2。水壩底部固定，頂部受到水的壓力。我們將使用平面應(yīng)變單元來(lái)分析水壩的變形和應(yīng)力分布。#導(dǎo)入必要的庫(kù)

importnumpyasnp

fromfenicsimport*

#定義幾何參數(shù)

height=50.0

width=10.0

E=30e9#彈性模量

nu=0.2#泊松比

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

mesh=RectangleMesh(Point(0,0),Point(width,height),10,50)

#定義邊界條件

defbottom_boundary(x,on_boundary):

returnnear(x[1],0.0)

deftop_boundary(x,on_boundary):

returnnear(x[1],height)

#定義材料屬性

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

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

#定義變分形式

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

du=TrialFunction(V)

u_=TestFunction(V)

u=Function(V)

#應(yīng)力張量

defsigma(v):

returnlmbda*tr(eps(v))*Identity(2)+2*mu*eps(v)

#應(yīng)變張量

defeps(v):

returnsym(nabla_grad(v))

#定義邊界條件

bc_bottom=DirichletBC(V,Constant((0,0)),bottom_boundary)

bc_top=DirichletBC(V.sub(1),Constant(-100000),top_boundary)#水的壓力

#定義外力

f=Constant((0,0))

#定義變分問(wèn)題

F=inner(sigma(du),eps(u_))*dx-inner(f,u_)*dx

#求解

solve(F==0,u,[bc_bottom,bc_top])

#輸出結(jié)果

file=File("displacement.pvd")

file<<u3.2軸對(duì)稱問(wèn)題3.2.1軸對(duì)稱問(wèn)題的定義軸對(duì)稱問(wèn)題是指結(jié)構(gòu)關(guān)于某一軸對(duì)稱，且所有物理量（如應(yīng)力、應(yīng)變和位移）都與該軸的徑向距離有關(guān)，而與角度無(wú)關(guān)。這類問(wèn)題常見(jiàn)于管道、圓柱體和圓盤等結(jié)構(gòu)的分析中。3.2.2軸對(duì)稱問(wèn)題的有限元模型在有限元分析中，軸對(duì)稱問(wèn)題可以通過(guò)建立二維模型來(lái)簡(jiǎn)化計(jì)算，其中模型沿徑向和軸向展開，而忽略角度方向的變化。這樣可以大大減少計(jì)算資源的需求，同時(shí)保持分析的準(zhǔn)確性。3.2.2.1示例：軸對(duì)稱問(wèn)題的有限元分析假設(shè)我們有一個(gè)圓柱形管道，其內(nèi)徑為50mm，外徑為100mm，長(zhǎng)度為1000mm，材料為鋼，彈性模量為200GPa，泊松比為0.3。管道內(nèi)受到10MPa的壓力。我們將使用軸對(duì)稱單元來(lái)分析管道的變形和應(yīng)力分布。#導(dǎo)入必要的庫(kù)

importnumpyasnp

fromfenicsimport*

#定義幾何參數(shù)

inner_radius=50.0

outer_radius=100.0

length=1000.0

E=200e9#彈性模量

nu=0.3#泊松比

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

mesh=Mesh()

editor=MeshEditor()

editor.open(mesh,"interval",2)

editor.init_vertices(100)

foriinrange(100):

editor.add_vertex(i,[i*(outer_radius-inner_radius)/99+inner_radius,0])

editor.close()

#定義邊界條件

definner_boundary(x,on_boundary):

returnnear(x[0],inner_radius)

defouter_boundary(x,on_boundary):

returnnear(x[0],outer_radius)

#定義材料屬性

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

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

#定義變分形式

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

du=TrialFunction(V)

u_=TestFunction(V)

u=Function(V)

#應(yīng)力張量

defsigma(v):

returnlmbda*tr(eps(v))*Identity(2)+2*mu*eps(v)

#應(yīng)變張量

defeps(v):

returnsym(nabla_grad(v))

#定義邊界條件

bc_inner=DirichletBC(V,Constant((0,0)),inner_boundary)

bc_outer=DirichletBC(V.sub(0),Constant(0),outer_boundary)

#定義外力

f=Constant((0,0))

p=Constant(-10e6)#內(nèi)部壓力

#定義變分問(wèn)題

F=inner(sigma(du),eps(u_))*dx-inner(f,u_)*dx-p*u_[0]*ds(1)

#求解

solve(F==0,u,[bc_inner,bc_outer])

#輸出結(jié)果

file=File("displacement.pvd")

file<<u3.3維彈性問(wèn)題3.3.1維問(wèn)題的單元類型三維彈性問(wèn)題需要使用三維單元來(lái)模擬，常見(jiàn)的三維單元包括四面體單元、六面體單元等。這些單元能夠準(zhǔn)確地描述三維空間中的應(yīng)力和應(yīng)變分布。3.3.2維問(wèn)題的有限元分析三維問(wèn)題的有限元分析比平面和軸對(duì)稱問(wèn)題復(fù)雜，因?yàn)樗枰紤]所有三個(gè)方向上的應(yīng)力和應(yīng)變。然而，現(xiàn)代有限元軟件能夠高效地處理這類問(wèn)題，使得分析三維結(jié)構(gòu)成為可能。3.3.2.1示例：三維問(wèn)題的有限元分析考慮一個(gè)立方體，其邊長(zhǎng)為100mm，材料為鋁，彈性模量為70GPa，泊松比為0.33。立方體的一側(cè)固定，另一側(cè)受到100N的力。我們將使用三維六面體單元來(lái)分析立方體的變形和應(yīng)力分布。#導(dǎo)入必要的庫(kù)

importnumpyasnp

fromfenicsimport*

#定義幾何參數(shù)

length=100.0

E=70e9#彈性模量

nu=0.33#泊松比

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

mesh=BoxMe