強(qiáng)度計(jì)算.數(shù)值計(jì)算方法：非線性分析：非線性動力學(xué)與振動

強(qiáng)度計(jì)算.數(shù)值計(jì)算方法：非線性分析：非線性動力學(xué)與振動1非線性動力學(xué)概述非線性動力學(xué)是研究非線性系統(tǒng)動力學(xué)行為的學(xué)科，它關(guān)注的是系統(tǒng)中力與位移、速度或加速度之間關(guān)系不遵循線性比例的物理現(xiàn)象。非線性動力學(xué)系統(tǒng)的行為往往復(fù)雜多變，包括但不限于混沌、分岔、周期性運(yùn)動和準(zhǔn)周期性運(yùn)動等。1.1非線性方程非線性系統(tǒng)的數(shù)學(xué)描述通常涉及非線性微分方程。例如，考慮一個簡單的非線性振動系統(tǒng)，如Duffing方程：Duffing方程:m*d^2x/dt^2+c*dx/dt+k*x+α*x^3=F(t)其中，m是質(zhì)量，c是阻尼系數(shù)，k是線性剛度系數(shù)，α是非線性剛度系數(shù)，F(xiàn)t1.2數(shù)值解法由于非線性方程的解析解往往難以找到，數(shù)值方法成為研究非線性動力學(xué)的主要工具。例如，使用Runge-Kutta方法求解Duffing方程：importnumpyasnp

fromegrateimportsolve_ivp

defduffing(t,y,m,c,k,alpha,F):

"""

定義Duffing方程的右端函數(shù)

y[0]=x,y[1]=dx/dt

"""

return[y[1],(F(t)-c*y[1]-k*y[0]-alpha*y[0]**3)/m]

#參數(shù)設(shè)置

m=1.0

c=0.1

k=1.0

alpha=-1.0

F=lambdat:np.cos(t)

#初始條件

y0=[0,0]

#時(shí)間范圍

t_span=(0,100)

#使用solve_ivp求解

sol=solve_ivp(duffing,t_span,y0,args=(m,c,k,alpha,F),dense_output=True)

#繪制結(jié)果

importmatplotlib.pyplotasplt

t=np.linspace(t_span[0],t_span[1],1000)

plt.plot(t,sol.sol(t)[0])

plt.xlabel('時(shí)間')

plt.ylabel('位移')

plt.title('Duffing方程的數(shù)值解')

plt.show()2振動理論基礎(chǔ)振動理論是研究物體在力的作用下產(chǎn)生周期性或非周期性運(yùn)動的學(xué)科。在非線性動力學(xué)中，振動理論被用來分析和預(yù)測非線性系統(tǒng)的動態(tài)響應(yīng)。2.1簡諧振動簡諧振動是最基本的振動形式，其運(yùn)動方程為：線性振動方程:m*d^2x/dt^2+c*dx/dt+k*x=02.2非簡諧振動非簡諧振動則涉及非線性項(xiàng)，如上述Duffing方程所示。非簡諧振動的分析通常需要更復(fù)雜的數(shù)學(xué)工具和數(shù)值方法。3非線性系統(tǒng)的分類非線性系統(tǒng)可以根據(jù)其非線性特性的不同進(jìn)行分類：硬彈簧系統(tǒng)：非線性剛度系數(shù)為正，系統(tǒng)在大位移時(shí)表現(xiàn)出更硬的特性。軟彈簧系統(tǒng)：非線性剛度系數(shù)為負(fù)，系統(tǒng)在大位移時(shí)表現(xiàn)出更軟的特性。阻尼非線性系統(tǒng)：阻尼力與速度的非線性關(guān)系，如平方律阻尼。外部激勵非線性系統(tǒng)：系統(tǒng)受到非線性外部力的作用，如沖擊力或隨機(jī)力。4非線性振動的特性非線性振動系統(tǒng)具有以下特性：多解性：對于相同的初始條件，非線性系統(tǒng)可能有多個穩(wěn)定解?；煦缧袨椋涸谀承l件下，非線性系統(tǒng)可能表現(xiàn)出混沌行為，即對初始條件的極端敏感性。分岔現(xiàn)象：隨著系統(tǒng)參數(shù)的變化，非線性系統(tǒng)的動態(tài)行為可能突然改變，這種現(xiàn)象稱為分岔。周期與準(zhǔn)周期運(yùn)動：非線性系統(tǒng)可能表現(xiàn)出復(fù)雜的周期性或準(zhǔn)周期性運(yùn)動，這些運(yùn)動可能無法通過簡單的線性分析預(yù)測。4.1分岔圖示例分岔圖是展示系統(tǒng)參數(shù)變化時(shí)，系統(tǒng)穩(wěn)定解的變化情況的圖形。以下是一個使用Python繪制Duffing方程分岔圖的例子：defduffing_map(y,t,m,c,k,alpha,F):

"""

定義Duffing方程的離散映射

"""

dt=0.01

returny+dt*[y[1],(F(t)-c*y[1]-k*y[0]-alpha*y[0]**3)/m]

#參數(shù)設(shè)置

m=1.0

c=0.1

k=1.0

alpha=np.linspace(-1.0,1.0,1000)

F=lambdat:np.cos(t)

#初始條件

y0=[0,0]

#時(shí)間范圍

t_span=(0,1000)

#繪制分岔圖

plt.figure(figsize=(10,5))

forainalpha:

sol=solve_ivp(duffing,t_span,y0,args=(m,c,k,a,F),dense_output=True)

t=np.linspace(t_span[0],t_span[1],1000)

plt.plot([a]*len(t),sol.sol(t)[0],'.',markersize=0.5)

plt.xlabel('非線性剛度系數(shù)')

plt.ylabel('位移')

plt.title('Duffing方程的分岔圖')

plt.show()以上內(nèi)容涵蓋了非線性動力學(xué)與振動的基礎(chǔ)理論、振動理論基礎(chǔ)、非線性系統(tǒng)的分類以及非線性振動的特性，包括具體的數(shù)值計(jì)算方法和代碼示例。5數(shù)值計(jì)算方法：非線性動力學(xué)與振動5.1有限元法簡介有限元法（FiniteElementMethod,FEM）是一種廣泛應(yīng)用于工程分析的數(shù)值方法，尤其在解決復(fù)雜的非線性動力學(xué)與振動問題中表現(xiàn)出色。它將連續(xù)的結(jié)構(gòu)或系統(tǒng)離散化為有限數(shù)量的單元，每個單元用一組節(jié)點(diǎn)來表示，通過在這些節(jié)點(diǎn)上求解微分方程的近似解，進(jìn)而得到整個結(jié)構(gòu)或系統(tǒng)的解。5.1.1原理有限元法基于變分原理和加權(quán)殘值法。在非線性動力學(xué)分析中，結(jié)構(gòu)的響應(yīng)可能受到材料非線性、幾何非線性或接觸非線性的影響。這些非線性效應(yīng)使得問題的解析解難以獲得，而有限元法通過將問題轉(zhuǎn)化為一系列線性或非線性代數(shù)方程組，可以有效地求解。5.1.2內(nèi)容離散化：將連續(xù)體劃分為有限個單元，每個單元用多項(xiàng)式函數(shù)近似。形函數(shù)：描述單元內(nèi)位移與節(jié)點(diǎn)位移之間的關(guān)系。剛度矩陣：基于單元的形函數(shù)和材料屬性，計(jì)算單元的剛度矩陣。非線性方程組：在考慮非線性效應(yīng)后，形成非線性方程組。求解算法：使用迭代方法求解非線性方程組，如Newton-Raphson法。5.2非線性方程的數(shù)值解法非線性方程的數(shù)值解法是解決非線性動力學(xué)與振動問題的關(guān)鍵。在有限元分析中，非線性方程組通常通過迭代方法求解。5.2.1原理迭代方法基于逐步逼近的策略，通過一系列的迭代步驟，逐步修正解的估計(jì)值，直到滿足收斂準(zhǔn)則。在非線性動力學(xué)分析中，常用的迭代方法有Newton-Raphson法、Arc-Length法等。5.2.2內(nèi)容Newton-Raphson法：基于泰勒級數(shù)展開，通過求解線性化后的方程組來逐步逼近非線性方程的解。Arc-Length法：在求解過程中引入一個控制參數(shù)，以確保在非線性路徑上穩(wěn)定地前進(jìn)。5.2.3示例代碼#Newton-Raphson法求解非線性方程

importnumpyasnp

deff(x):

#定義非線性方程

returnx**3-2*x-5

defdf(x):

#定義非線性方程的導(dǎo)數(shù)

return3*x**2-2

defnewton_raphson(x0,tol=1e-6,max_iter=100):

#Newton-Raphson迭代求解

x=x0

foriinrange(max_iter):

x_new=x-f(x)/df(x)

ifabs(x_new-x)<tol:

returnx_new

x=x_new

returnNone

#初始猜測值

x0=2.0

#迭代求解

solution=newton_raphson(x0)

print("Solution:",solution)5.3時(shí)間積分方法時(shí)間積分方法用于求解隨時(shí)間變化的非線性動力學(xué)問題，如結(jié)構(gòu)的動力響應(yīng)分析。5.3.1原理時(shí)間積分方法通過在時(shí)間域內(nèi)離散化動力學(xué)方程，將連續(xù)的時(shí)間問題轉(zhuǎn)化為一系列離散的時(shí)間步問題。常見的方法有顯式和隱式積分方法，如Euler法、Runge-Kutta法、Newmark法等。5.3.2內(nèi)容顯式積分方法：如Euler法，計(jì)算速度快，但穩(wěn)定性條件嚴(yán)格。隱式積分方法：如Newmark法，穩(wěn)定性好，但計(jì)算量大。5.3.3示例代碼#使用Euler法進(jìn)行時(shí)間積分

importnumpyasnp

defacceleration(t,x,v):

#定義加速度函數(shù)

return-x+np.sin(t)

defeuler_integration(x0,v0,t0,tf,dt):

#Euler積分求解

t=t0

x=x0

v=v0

whilet<tf:

a=acceleration(t,x,v)

v=v+a*dt

x=x+v*dt

t=t+dt

returnx,v

#初始條件

x0=0.0

v0=0.0

t0=0.0

tf=10.0

dt=0.1

#時(shí)間積分求解

x_final,v_final=euler_integration(x0,v0,t0,tf,dt)

print("Finalposition:",x_final)

print("Finalvelocity:",v_final)5.4非線性動力學(xué)的數(shù)值模擬非線性動力學(xué)的數(shù)值模擬結(jié)合了有限元法、非線性方程的數(shù)值解法和時(shí)間積分方法，用于預(yù)測和分析復(fù)雜結(jié)構(gòu)在非線性動力學(xué)載荷下的響應(yīng)。5.4.1原理通過有限元法建立結(jié)構(gòu)的非線性動力學(xué)模型，然后使用非線性方程的數(shù)值解法和時(shí)間積分方法求解模型，得到結(jié)構(gòu)的動力響應(yīng)。5.4.2內(nèi)容模型建立：使用有限元法建立結(jié)構(gòu)的非線性動力學(xué)模型。求解：結(jié)合非線性方程的數(shù)值解法和時(shí)間積分方法求解模型。結(jié)果分析：分析結(jié)構(gòu)的動力響應(yīng)，如位移、速度、加速度等。5.4.3示例代碼#使用有限元法和時(shí)間積分方法進(jìn)行非線性動力學(xué)模擬

importnumpyasnp

deffem_nonlinear_dynamics(model,load,dt,tf):

#定義非線性動力學(xué)模擬函數(shù)

t=0.0

whilet<tf:

#更新載荷

model.update_load(load(t))

#求解非線性方程組

solution=model.solve_nonlinear()

#時(shí)間積分

model.time_integration(dt)

t+=dt

returnmodel.response()

#創(chuàng)建有限元模型

model=FiniteElementModel()

#定義載荷函數(shù)

defload(t):

returnnp.sin(t)

#時(shí)間步和總時(shí)間

dt=0.1

tf=10.0

#進(jìn)行非線性動力學(xué)模擬

response=fem_nonlinear_dynamics(model,load,dt,tf)

print("Response:",response)請注意，上述代碼示例是簡化版的，實(shí)際的有限元分析和非線性動力學(xué)模擬會涉及更復(fù)雜的數(shù)學(xué)模型和算法實(shí)現(xiàn)。6非線性分析6.1材料非線性分析6.1.1原理材料非線性分析關(guān)注材料在大應(yīng)變、大應(yīng)力或溫度變化下的行為。材料的非線性特性通常包括塑性、粘彈性、超彈性、蠕變和熱力學(xué)非線性。在塑性分析中，材料在達(dá)到屈服點(diǎn)后會發(fā)生永久變形，這需要使用塑性理論來描述。粘彈性材料在應(yīng)力作用下會隨時(shí)間發(fā)生變形，其行為可以用Kelvin-Voigt模型或Maxwell模型來模擬。超彈性材料，如橡膠和生物材料，具有高度非線性的應(yīng)力-應(yīng)變關(guān)系，通常使用Mooney-Rivlin或Ogden模型來描述。蠕變分析涉及材料在長時(shí)間應(yīng)力作用下的變形，而熱力學(xué)非線性則考慮溫度變化對材料性能的影響。6.1.2內(nèi)容材料非線性分析的核心是建立準(zhǔn)確的材料模型，并將其應(yīng)用于有限元分析中。例如，塑性分析中常用的模型有vonMises屈服準(zhǔn)則和Tresca屈服準(zhǔn)則。在有限元軟件中，如ABAQUS，可以定義材料屬性，包括彈性模量、泊松比、屈服強(qiáng)度和硬化參數(shù)，以進(jìn)行塑性分析。示例在ABAQUS中定義一個簡單的塑性材料模型：#ABAQUSPythonScriptfordefiningaplasticmaterial

fromabaqusimport*

fromabaqusConstantsimport*

fromodbAccessimport*

frommaterialimport*

fromsectionimport*

session.journalOptions.setValues(replayGeometry=COORDINATE,recoverGeometry=COORDINATE)

#Createanewmaterial

myModel=mdb.models['Model-1']

myMaterial=myModel.Material(name='PlasticMaterial')

#Defineelasticproperties

myMaterial.Elastic(table=((200e9,0.3),))

#Defineplasticproperties

myMaterial.Plastic(table=((250e6,0.0),(300e6,0.002),(350e6,0.005)))

#Assignmaterialtoasection

mySection=myModel.HomogeneousSolidSection(name='PlasticSection',material='PlasticMaterial',thickness=None)這段代碼定義了一個具有塑性特性的材料，并將其分配給一個固體截面。材料的彈性模量為200GPa，泊松比為0.3，塑性行為由屈服強(qiáng)度和硬化參數(shù)的表格定義。6.2幾何非線性分析6.2.1原理幾何非線性分析考慮結(jié)構(gòu)在大變形下的行為，其中變形會導(dǎo)致結(jié)構(gòu)幾何形狀的顯著變化，從而影響其力學(xué)性能。這種分析通常在結(jié)構(gòu)承受大位移、大旋轉(zhuǎn)或大應(yīng)變時(shí)進(jìn)行，例如在碰撞、成型或結(jié)構(gòu)穩(wěn)定性分析中。幾何非線性分析的關(guān)鍵是使用非線性方程求解器，如Newton-Raphson方法，來迭代求解結(jié)構(gòu)的平衡狀態(tài)。6.2.2內(nèi)容在幾何非線性分析中，結(jié)構(gòu)的變形狀態(tài)需要在每個時(shí)間步或載荷步進(jìn)行更新。這涉及到使用更新的位移和旋轉(zhuǎn)來重新計(jì)算結(jié)構(gòu)的幾何形狀，然后基于新的幾何形狀求解結(jié)構(gòu)的平衡狀態(tài)。在有限元分析中，這通常通過選擇適當(dāng)?shù)姆治鲱愋秃颓蠼馄髟O(shè)置來實(shí)現(xiàn)。示例在ABAQUS中設(shè)置幾何非線性分析：#ABAQUSPythonScriptforsettingupgeometricnonlinearity

fromabaqusimport*

fromabaqusConstantsimport*

fromstepimport*

session.journalOptions.setValues(replayGeometry=COORDINATE,recoverGeometry=COORDINATE)

#Createanewstepforgeometricnonlinearity

myStep=myModel.StaticStep(name='GeometricNonlinearityStep',previous='Initial',initialInc=0.1,maxNumInc=1000,nlgeom=ON)

#Settheanalysistypetogeometricnonlinearity

myStep.setValues(nlgeom=ON)這段代碼創(chuàng)建了一個新的靜力分析步驟，并設(shè)置了非線性幾何選項(xiàng)。nlgeom=ON表示分析將考慮幾何非線性。6.3接觸非線性分析6.3.1原理接觸非線性分析研究兩個或多個物體在接觸時(shí)的相互作用。接觸可以是滑動、粘著或分離的，這取決于接觸面的法向和切向行為。接觸分析通常需要定義接觸對，包括主面和從面，以及接觸屬性，如摩擦系數(shù)和接觸壓力。6.3.2內(nèi)容接觸非線性分析在工程設(shè)計(jì)中非常重要，尤其是在摩擦、磨損和碰撞問題中。在有限元分析中，接觸分析通常需要定義接觸對，設(shè)置接觸屬性，并選擇適當(dāng)?shù)慕佑|算法，如罰函數(shù)法或拉格朗日乘子法。示例在ABAQUS中定義接觸對：#ABAQUSPythonScriptfordefiningacontactpair

fromabaqusimport*

fromabaqusConstantsimport*

frominteractionimport*

session.journalOptions.setValues(replayGeometry=COORDINATE,recoverGeometry=COORDINATE)

#Createanewinteractionforcontact

myContact=myModel.ContactExp(name='ContactInteraction')

#Definethecontactpair

myContact.ContactProperty('ContactProperty')

myContact.TangentialBehavior(formulation=FINITE,directionality=ISOTROPIC,slipRateDependency=OFF,pressureDependency=OFF,temperatureDependency=OFF,dependencies=0,table=((0.3,),))

myContact.NormalBehavior(pressureOverclosure=HARD,allowSeparation=ON,generateContactElements=ON)

#Assignthecontactpair

myContact=myMeractions['ContactInteraction']

myContact.includedPairs.setValuesInStep(stepName='Initial',useAllstar=ON)這段代碼定義了一個接觸對，并設(shè)置了接觸屬性，包括摩擦系數(shù)和接觸行為。6.4多物理場耦合非線性分析6.4.1原理多物理場耦合非線性分析考慮了不同物理場之間的相互作用，如結(jié)構(gòu)力學(xué)、熱力學(xué)、電磁學(xué)和流體力學(xué)。在這些分析中，一個物理場的變化會影響另一個物理場，反之亦然。例如，在熱-結(jié)構(gòu)耦合分析中，溫度變化會導(dǎo)致材料性能的變化，從而影響結(jié)構(gòu)的力學(xué)行為；同時(shí)，結(jié)構(gòu)的變形也會影響熱傳導(dǎo)路徑，進(jìn)而影響溫度分布。6.4.2內(nèi)容多物理場耦合非線性分析通常需要在有限元軟件中定義多個物理場，并設(shè)置它們之間的耦合關(guān)系。這可能涉及到定義材料的多物理場屬性，如熱膨脹系數(shù)、熱導(dǎo)率和電磁性能，以及設(shè)置載荷和邊界條件，以反映不同物理場的相互作用。示例在ABAQUS中設(shè)置熱-結(jié)構(gòu)耦合分析：#ABAQUSPythonScriptforsettingupthermo-mechanicalcoupling

fromabaqusimport*

fromabaqusConstantsimport*

fromstepimport*

fromloadimport*

fromboundaryConditionimport*

session.journalOptions.setValues(replayGeometry=COORDINATE,recoverGeometry=COORDINATE)

#Createanewstepforthermo-mechanicalcoupling

myStep=myModel.CoupledTempDisplacementStep(name='ThermoMechanicalStep',previous='Initial',initialInc=0.1,maxNumInc=1000,nlgeom=ON)

#Definethermalload

myModel.ConcentratedHeatFlux(name='HeatLoad',createStepName='ThermoMechanicalStep',region=myModel.rootAssembly.sets['HeatLoadRegion'],cf1=1000.0)

#Definemechanicalboundarycondition

myModel.DisplacementBC(name='MechanicalBC',createStepName='ThermoMechanicalStep',region=myModel.rootAssembly.sets['MechanicalBCRegion'],u1=0.0,u2=0.0,u3=0.0,ur1=UNSET,ur2=UNSET,ur3=UNSET,amplitude=UNSET,fixed=OFF,distributionType=UNIFORM,fieldName='',localCsys=None)這段代碼創(chuàng)建了一個熱-結(jié)構(gòu)耦合分析步驟，并定義了熱載荷和機(jī)械邊界條件。熱載荷施加在名為HeatLoadRegion的區(qū)域，而機(jī)械邊界條件則施加在MechanicalBCRegion區(qū)域，限制了所有方向的位移。以上示例展示了如何在ABAQUS中進(jìn)行材料非線性、幾何非線性、接觸非線性和多物理場耦合非線性分析的設(shè)置。這些分析對于理解和預(yù)測復(fù)雜工程結(jié)構(gòu)在各種載荷條件下的行為至關(guān)重要。7非線性動力學(xué)與振動應(yīng)用7.1非線性振動的工程實(shí)例在工程領(lǐng)域，非線性振動現(xiàn)象普遍存在，尤其是在大型結(jié)構(gòu)、橋梁、風(fēng)力發(fā)電機(jī)葉片等復(fù)雜系統(tǒng)中。這些系統(tǒng)由于材料的非線性、幾何非線性或外部激勵的非線性，其振動特性不能用線性理論準(zhǔn)確描述。例如，風(fēng)力發(fā)電機(jī)葉片在高速旋轉(zhuǎn)時(shí)，由于離心力的作用，其剛度會發(fā)生變化，這種現(xiàn)象就屬于幾何非線性。7.1.1示例：風(fēng)力發(fā)電機(jī)葉片的非線性振動分析假設(shè)我們有一臺風(fēng)力發(fā)電機(jī)葉片，其長度為L，質(zhì)量為m，剛度為k，且剛度隨離心力的變化而變化，可以表示為k=k0+αω^2，其中k0是靜止時(shí)的剛度，α是與葉片材料和幾何形狀相關(guān)的常數(shù)，ω是葉片的旋轉(zhuǎn)角速度。我們可以通過數(shù)值方法，如有限元法，來分析葉片在不同旋轉(zhuǎn)速度下的振動特性。importnumpyasnp

fromegrateimportsolve_ivp

#定義參數(shù)

m=1000#質(zhì)量，單位：kg

k0=1e6#初始剛度，單位：N/m

alpha=1e3#剛度變化系數(shù)，單位：N/(m*rad^2)

omega=10#旋轉(zhuǎn)角速度，單位：rad/s

#定義非線性振動方程

defnonlinear_vibration(t,y):

k=k0+alpha*omega**2#計(jì)算剛度

dydt=[y[1],-(k/m)*y[0]]#一階微分方程組

returndydt

#初始條件

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

#時(shí)間范圍

t_span=(0,10)

#解微分方程

sol=solve_ivp(nonlinear_vibration,t_span,y0,t_eval=np.linspace(0,10,100))

#輸出結(jié)果

print(sol.t)#時(shí)間點(diǎn)

print(sol.y[0])#位移7.2非線性動力學(xué)在結(jié)構(gòu)工程中的應(yīng)用非線性動力學(xué)在結(jié)構(gòu)工程中的應(yīng)用廣泛，包括地震工程、橋梁工程、高層建筑的風(fēng)振分析等。在這些應(yīng)用中，結(jié)構(gòu)的非線性響應(yīng)是設(shè)計(jì)和評估安全性的關(guān)鍵因素。7.2.1示例：地震作用下橋梁的非線性動力學(xué)分析考慮一座橋梁在地震作用下的非線性動力學(xué)響應(yīng)。橋梁的簡化模型可以視為一個帶有非線性彈簧的單自由度系統(tǒng)。地震加速度可以通過一個時(shí)間序列來模擬，我們使用Newmark-β方法來求解橋梁的響應(yīng)。importnumpyasnp

#定義參數(shù)

m=1e3#質(zhì)量，單位：kg

k0=1e7#初始剛度，單位：N/m

c=1e4#阻尼系數(shù)，單位：N/(m/s)

alpha=1e3#剛度變化系數(shù)，單位：N/(m*rad^2)

omega=1#地震頻率，單位：rad/s

A=1#地震加速度幅值，單位：m/s^2

t_end=10#分析時(shí)間，單位：s

dt=0.01#時(shí)間步長，單位：s

#地震加速度時(shí)間序列

t=np.arange(0,t_end,dt)

a=A*np.sin(omega*t)

#Newmark-β方法參數(shù)

beta=0.25

gamma=0.5

#初始化

u=np.zeros_like(t)

v=np.zeros_like(t)

a_sys=np.zeros_like(t)

u[0]=0.1#初始位移

v[0]=0#初始速度

#Newmark-β方法求解

foriinrange(1,len(t)):

t_n=t[i-1]

t_np1=t[i]

u_n=u[i-1]

v_n=v[i-1]

a_n=a[i-1]

#預(yù)測位移和速度

u_np1_pred