用Python进行线性和非线性结构分析

01

案例介绍

这个案例中主要介绍CivilFEM2018高级功能——Python脚步在结构建模、加载、分析与结构后处理各个过程中的应用。

案例分析了由杆件和弹簧组成的简单结构，杆件一端固定，另一端与一线性弹簧连接，这端可以在垂直方向上移动，并受到30磅的集中荷载P，方向如下图，加载随时间变化而变化，最后达到30磅。

02

结构基本参数

杆件截面：实心圆柱

杆件弹性模量： 1e7

杆件泊松比： 0.3

线性弹簧弹性系数： 6lb/in

03

使用Python构建结构模型

使用Python建模、定义截面属性、材质、结构元素和划分网格。

newDocument("Dim2D","Static","Structural")

ConfigUnits.System = "BInch"

import matplotlib.pyplot as pl

import matplotlib.animation as animation

# Geometry

createPoint("1",Point2d(0,0))

createPoint("2",Point2d(100,1))

createLine("Truss","1","2")

# Materials and sections

createGenericMat("Truss_mat",Double(1e+07),Double(0.3),Double(0),Double(0))

createGenericSection("Truss_sect","Truss_mat",Double(1),Double(1e-10),Double(1e-10),Double(0))

# Structural elements

createSETruss("Truss","Truss_sect","Truss_sect","Truss")

structuralElem0 = StructuralElementsContainer.Find("Truss")

structuralElem0.MeshTool1D.ParameterMesh.NumberOfSegments = 1

createPointSpring("Spring","ToGround",Double(6),"Truss",Point2d(100,1),False,True)

makeZoom("Fit")

# Mesh

mesh()

04

施加荷载、定义边界条件和荷载工况

施加时序荷载：

1s: 5lb; 4s: 8lb; 4.8s: 12 lb; 6.6s: 30lb

定义杆件两端约束；

创建4个荷载工况。

# Loads and BCs

createLoadGroup("Load_1")

loadGroups4 = LoadGroupsContainer.Find("Load_1")

addPointLoad([loadGroups4],"Load","Truss",Point2d(100,1),Vector2d(0,-1),Double(5))

createLoadGroup("Load_4")

loadGroups4 = LoadGroupsContainer.Find("Load_4")

addPointLoad([loadGroups4],"Load","Truss",Point2d(100,1),Vector2d(0,-1),Double(8))

createLoadGroup("Load_4_8")

loadGroups4 = LoadGroupsContainer.Find("Load_4_8")

addPointLoad([loadGroups4],"Load","Truss",Point2d(100,1),Vector2d(0,-1),Double(12))

createLoadGroup("Load_6_6")

loadGroups4 = LoadGroupsContainer.Find("Load_6_6")

addPointLoad([loadGroups4],"Load","Truss",Point2d(100,1),Vector2d(0,-1),Double(30))

createBCGroup("BC")

bc5 = BCGroupsContainer.Find("BC")

addPointBC([bc5],"1","Truss",Point2d(0,0),True,True)

addPointBC([bc5],"2","Truss",Point2d(100,1),True,False)

# LCs

createLoadCase("LC_1_0s")

loadCase1 = LoadCasesContainer.Find("LC_1_0s")

addLoadGroupToLoadCase([loadCase1],"Load_1",Double(1))

addBCGroupToLoadCase([loadCase1],"BC")

createLoadCase("LC_4_0s")

loadCase1 = LoadCasesContainer.Find("LC_4_0s")

addLoadGroupToLoadCase([loadCase1],"Load_4",Double(1))

addBCGroupToLoadCase([loadCase1],"BC")

loadCase1.LoadCaseCalculationTime = Double(4.0)

createLoadCase("LC_4_8s")

loadCase1 = LoadCasesContainer.Find("LC_4_8s")

addLoadGroupToLoadCase([loadCase1],"Load_4_8",Double(1))

addBCGroupToLoadCase([loadCase1],"BC")

loadCase1.LoadCaseCalculationTime = Double(4.8)

createLoadCase("LC_6_6s")

loadCase1 = LoadCasesContainer.Find("LC_6_6s")

addLoadGroupToLoadCase([loadCase1],"Load_6_6",Double(1))

addBCGroupToLoadCase([loadCase1],"BC")

loadCase1.LoadCaseCalculationTime = Double(6.6)

05

线性求解分析

solve()

# Obtain nodal results for linear case and store them in a list

D_lin=[0.0]*5

Node=Double(0.0)

linear=getResults(r"LC_1_0s.rcf")

linear.getNodeResult(Node, "UTy", 2)

D_lin[1]=-(Node.getValue())

linear=getResults(r"LC_4_0s.rcf")

linear.getNodeResult(Node, "UTy", 2)

D_lin[2]=-(Node.getValue())

linear=getResults(r"LC_4_8s.rcf")

linear.getNodeResult(Node, "UTy", 2)

D_lin[3]=-(Node.getValue())

linear=getResults(r"LC_6_6s.rcf")

linear.getNodeResult(Node, "UTy", 2)

D_lin[4]=-(Node.getValue())

06

非线性求解分析

非线性求解主要聚焦结构的大位移，在Python脚本中对此结构出现的大位移进行分析求解。

# Add individual solution controls to the LCs for the non linear case

loadCase1 = LoadCasesContainer.Find("LC_1_0s")

loadCase1.SolutionControl.IsParticularControl = True

loadCase1.SolutionControl.InitialTimeStep = Formula("0.1")

loadCase1.SolutionControl.MinimumTimeStep = Formula("0.1")

loadCase1.SolutionControl.MaximumTimeStep = Formula("0.1")

loadCase1 = LoadCasesContainer.Find("LC_4_0s")

loadCase1.SolutionControl.IsParticularControl = True

loadCase1.SolutionControl.InitialTimeStep = Formula("0.1/3")

loadCase1.SolutionControl.MinimumTimeStep = Formula("0.1/3")

loadCase1.SolutionControl.MaximumTimeStep = Formula("0.1/3")

loadCase1 = LoadCasesContainer.Find("LC_4_8s")

loadCase1.SolutionControl.IsParticularControl = True

loadCase1.SolutionControl.InitialTimeStep = Formula("0.125")

loadCase1.SolutionControl.MinimumTimeStep = Formula("0.125")

loadCase1.SolutionControl.MaximumTimeStep = Formula("0.125")

loadCase1 = LoadCasesContainer.Find("LC_6_6s")

loadCase1.SolutionControl.IsParticularControl = True

loadCase1.SolutionControl.InitialTimeStep = Formula("1/18")

loadCase1.SolutionControl.MinimumTimeStep = Formula("1/18")

loadCase1.SolutionControl.MaximumTimeStep = Formula("1/18")

ConfigSolver.CreateIntermediateNonLinearResults = True

ConfigSolver.LargeDeflections = True

solve()

# Add non-linear results to vector

D_n_lin=[0.0]*67

for i in range(1,11):

non_linear=getResults(r"LC_1_0s_"+str(i)+".rcf")

non_linear.getNodeResult(Node, "UTy", 2)

D_n_lin[i]=-(Node.getValue())

for i in range(11,41):

non_linear=getResults(r"LC_4_0s_"+str(i)+".rcf")

non_linear.getNodeResult(Node, "UTy", 2)

D_n_lin[i]=-(Node.getValue())

for i in range(41,49):

non_linear=getResults(r"LC_4_8s_"+str(i)+".rcf")

non_linear.getNodeResult(Node, "UTy", 2)

D_n_lin[i]=-(Node.getValue())

for i in range(49,67):

non_linear=getResults(r"LC_6_6s_"+str(i)+".rcf")

non_linear.getNodeResult(Node, "UTy", 2)

D_n_lin[i]=-(Node.getValue())

07

结果后处理

在图形界面中比较节点2在线性和非线性分析两种情况下位移结果。

# Plot results

fig=pl.figure(figsize=(30, 10), dpi=72)

ax1 = pl.subplot(1, 3, 1)

F_lin=[0, 5, 8, 12, 30]

ax1.plot(D_lin, F_lin, marker='o', color="blue", linewidth=1.0, linestyle="-", label='Linear')

F_n_lin=[]

for x in range (0,11):

F_n_lin.append(5.0/10.0*x)

for x in range (1,31):

F_n_lin.append(3.0/30.0*x+5.0)

for x in range (1,9):

F_n_lin.append(4.0/8.0*x+8.0)

for x in range (1,19):

F_n_lin.append(18.0/18.0*x+12.0)

# First subplot

ax1.plot(D_n_lin, F_n_lin, marker='.', color="red", linewidth=1.0, linestyle="-", label='Non-linear')

ax1.legend(loc='upper left')

ax1.set_ylabel('Force(lbf)')

ax1.set_xlabel('Abs displacement (in)')

ax1.set_title('Linear VS Non-linear')

ax1.set_ylim(0.0, 30.3)

ax1.set_xlim(0.0, 2.7)

ax1.grid(True)

动画模拟杆件在非线性分析工况中随着时序加载过程的大位移变化过程。

# Animation of second subplot.

ax2 = pl.subplot(1, 3, 2)

ax2.plot([0,100], [0,1], marker='.', color="blue", linewidth=1.0, linestyle="-", label='No Load')

line, =ax2.plot([0,100], [0,1-D_n_lin[i]], marker='.', color="red", linewidth=1.0, linestyle="-", label='Non-linear deformation')

def animate(i):

line.set_data([0,100], [0,1-D_n_lin[i]])

return line,

ax2.set_title('Truss displacement evolution, non-linear case')

ax2.set_ylim(-1.8, 1.1)

ax2.set_xlim(0.0, 101)

ax2.set_ylabel('Real position (in)')

ax2.set_xlabel('Length of truss (in)')

ax2.grid(True)

line_ani = animation.FuncAnimation(fig, animate, frames=67, interval=5, blit=True)

可见，CivilFEM2018整个分析过程都可以使用Python脚本，对于熟悉Python的用户，这个功能可以大幅提升工作效率。

CivilFEM 2018 powered by Marc （简称CivilFEM）是一款高级非线性土木工程专用有限元商业软件，由西班牙INGECIBER公司开发，该公司同时开发ANSYS/CIVILFEM土木模块。独立的CivilFEM商业软件采用了Marc非线性求解器，方便易用的建模环境，快速高效的求解过程，丰富的后处理功能，CivilFEM2018正在成为土木工程师必备的卓越软件工具。上海图杰信息科技公司是CivilFEM2018在中国地区的销售和技术支持机构。