基于Python的结构力学位移法编程求解

最新推荐文章于 2024-04-29 11:00:38 发布

qq_45885422

最新推荐文章于 2024-04-29 11:00:38 发布

阅读量3.6k

点赞数 2

文章标签： python

本文链接：https://blog.csdn.net/qq_45885422/article/details/121977389

版权

如果想用Fortran编，可以参考基于Fortran的结构力学位移法编程求解

1.背景介绍

位移法的基本思路使以结点位移(广义位移)作为基本未知量，写出由未知位移表示的应变，由物理方程写出仍由未知位移表示的应力表达式，最后由平衡条件解出所有的未知位移。

在计算机科学飞速发展的今天，适合于计算机应用的“有限元法”正在逐步取代其他方法而成为结构分析方法的主流，并已发展为一门独立的新兴学科。

本代码可用于带符号运算的平面桁架结构的计算。

2.编程思路

2.1 结点编号、元素编号，建立总体坐标系xoy

此步一般由题目直接给出，本步骤以下图为例。

已知条件为

$u_{1}=v_{1}=u_{4}=v_{4}=0$

$u_{3}=b$

$P_{2x}=X_{2},P_{2y}=Y_{2},P_{3x}=0,P_{3y}=0$

桁架各单元的尺寸和倾角列于下表。

2.2 建立各元素在总体坐标系xoy 中的刚度矩阵

记 $\lambda =cos\theta ,\mu =sin\theta$

则杆单元为ij的刚度矩阵为：

2.3 建立总刚度方程

当全结构有N个结点，将N个结点的平衡方程联立起来，组成全结构的平衡方程为：

首先按编号顺序列出结点位移和结点力列阵：

$\left \{ \delta \right \}=\begin{bmatrix} u_{1} &v_{1} &u_{2} &v_{2} & u_{3} & v_{3} & u_{4} & v_{4} \end{bmatrix}'$

$\left \{ P\right \}=\begin{bmatrix} P_{1x} &P_{1y} &P_{2x} &P_{2y}&P_{3x} &P_{3y} &P_{4x} &P_{4y} \end{bmatrix}'$

接着组集总刚度矩阵：

$K_{ij}=\frac{EA}{L_{ij}}\begin{bmatrix} \lambda ^{2} &\lambda \mu &-\lambda ^{2}&-\lambda \mu \\ \lambda \mu&\mu^{2} & -\lambda \mu &-\mu^{2} \\ -\lambda ^{2} &-\lambda \mu &\lambda ^{2}&\lambda \mu \\ -\lambda \mu&-\mu^{2} & \lambda \mu &\mu^{2} \end{bmatrix}=\begin{bmatrix} \begin{bmatrix} K_{ii} \end{bmatrix} &\begin{bmatrix} K_{ij} \end{bmatrix} \\ \begin{bmatrix} K_{ji} \end{bmatrix} & \begin{bmatrix} K_{jj} \end{bmatrix} \end{bmatrix}$

2.4 消除刚体位移

总刚度方程在没有引入位移约束条件之前，是不能求解的。因此在结构分析中，采用引入边界条件的办法删去相关的行(列)，之后的总刚度矩阵[ K ]变成可逆的，此时总刚度方程才能求解。

引入位移约束条件：

$u_{1}=v_{1}=u_{4}=v_{4}=0$

接着就可求解方程了，由计算机完成即可，略。

2.5 求解杆轴力

$N_{ij}=\frac{EA}{L_{ij}}\begin{bmatrix} -\lambda& -\mu&\lambda& \mu \end{bmatrix}\begin{Bmatrix} u_{i}\\ v_{i}\\ u_{j}\\ v_{j} \end{Bmatrix}$

3.代码介绍

3.1 输入输出部分

为了方便每次输入数据，因此使用txt文件作为输入输出文件，可以避免输出数据后多次重复输入的麻烦。

输入文件格式如下：

第1部分：杆件数结点数，如杆件数6，节点数为4，写做：6 4
第2部分：杆结点i 杆结点j 倾角长度，如长度为1的23号杆倾角为45度，写做：2 3 45 1
第三部分：结点横向位移结点纵向位移结点横向力结点纵向力，如某节点受力为（Px,Py），位移为（0，0），写做：0 0 Px Py

以之前的数据为例，该输入文件写入如下：

注意：
书写格式时需全在英文状态下
输入量中间用空格隔开
结点位移未知时可用任意非零数替代
仅支持结点力为未知量时的计算

在运行代码时，直接输入文件的文件名（带后缀）就可实现求解。

相应的，输出结果保存在“计算结果.txt”中，如下所示：

3.2 完整代码

from numpy import *
from math import *
import sys,os,sympy,codecs

def main():
    path=input('请输入文件名：')
    try:
        file=codecs.open(path,'r',encoding='utf-8')
    except:
        print('文件不存在，请检查，程序结束运行。\n')
    data=file.readlines()
    time=0
    line=data[time].strip('\r\n')
    temps=line.split(' ',line.count(' '))
    num=int(temps[0])
    point=int(temps[1])

    m=list(range(num))
    mm=zeros((point*2,point*2))
    uv=zeros(point*2)
    forces=[0]*(point*2)
    angle=zeros(num)
    length=zeros(num)
    num1=zeros(num)
    num2=zeros(num)

    time=time+1

    for i in range(num):
        line=data[time].strip('\r\n')
        temps=line.split(' ',line.count(' '))
        angle[i]=float(temps[2])
        length[i]=eval(temps[3])
        num1[i]=int(temps[0])
        num2[i]=int(temps[1])
        m[i]=mat(angle[i],length[i])
        # print(m[i],end='\n')#输出刚度矩阵
        m[i]=rebig(m[i],num1[i],num2[i],point)
        # print(m[i],end='\n')#输出扩阶后的刚度矩阵
        time=time+1

    for i in range(num):
        mm=mm+m[i]
    #print(mm,end='\n')#输出总刚度矩阵
    for i in range(point):
        line=data[time].strip('\r\n')
        temps=line.split(' ',line.count(' '))


        uv[i*2]=float(temps[0])
        uv[i*2+1]=float(temps[1])
        forces[i*2]=for_num(temps[2])
        forces[i*2+1]=for_num(temps[3])
        time=time+1
    file.close()        
    location=get_location_in_list(uv.tolist(),0)
    mms=delete(mm,[location], axis=1) 
    mms=delete(mms,[location], axis=0)
    if linalg.det(mms)<0.00000000001:
        print('输入数据出错',end='\n')
        os._exit(0)
    forces_delete=delete(forces,[location], axis=0)

    uv=add_add(point,uv,((mms.I).dot(forces_delete.T)).T,location)
    bar_force(angle,uv,length,num1,num2,num)
    uv=matrix(uv)
    forces=(mm.dot(uv.T))
    
    with open('计算结果.txt','a') as f:
        f.write('各结点力为：\n')
        for i in range(point):
            f.write('结点%d力：'%(i+1))
            f.write(str(forces[2*i,0]))
            f.write('  ,  ')
            f.write(str(forces[2*i+1,0]))
            f.write('\n')
        f.write('\n各结点位移为：\n')
        for i in range(point):
            f.write('结点%d位移：L/EA('%(i+1))
            f.write(str(uv[0,2*i],))
            f.write('  ,  ')
            f.write(str(uv[0,2*i+1]))
            f.write(')\n')
        f.close()
    #for i in range(point):
        #print('结点%d力：'%(i+1),forces[2*i,0],',',forces[2*i+1,0],end='\n')
    #for i in range(point):
        #print('结点%d位移：L/EA('%(i+1),uv[0,2*i],',',uv[0,2*i+1],end=')\n')

    print('计算完成，结果保存在"计算结果.txt"\n')


#建立元素刚度矩阵
def mat(angle,length):
    lamda=cos(angle/180*pi)
    miu=sin(angle/180*pi)
    return (matrix([[lamda**2,miu*lamda,-lamda**2,-lamda*miu],
                         [miu*lamda,miu**2,-lamda*miu,-miu**2],
                         [-lamda**2,-miu*lamda,lamda**2,lamda*miu],
                         [-lamda*miu,-miu**2,lamda*miu,miu**2]]))/length


#组装总刚度矩阵
def rebig(m,num1,num2,point):
    for i in range(point*2):
        if i not in [num1*2-2,num1*2-1,num2*2-2,num2*2-1]:
            m=insert(m,[i],zeros(m.shape[1]),axis=0)
    for i in range(point*2):
        if i not in [num1*2-2,num1*2-1,num2*2-2,num2*2-1]:      
            m=insert(m,[i],array([zeros(m.shape[0])]).T,axis=1)
    return m

def get_location_in_list(x, target):
    step = -1
    items = list()
    for i in range(x.count(target)):
        y = x[step + 1:].index(target)
        step = step + y + 1
        items.append(step)
    return items

def add_add(point,dis,add,loca):
    temp=[0]*(point*2)
    t=0
    for i in range(point*2):
         if i not in loca:
             temp[i]=add[t,0]
             t=t+1
         else:
             temp[i]=0
    return temp

def bar_force(angle,uv,length,num1,num2,num):
    with open('计算结果.txt','w') as f:
        f.write('各杆力为：\n')
        for i in range(num):
            lamda=cos(angle[i]/180*pi)
            miu=sin(angle[i]/180*pi)      
            bar_force=(-lamda*uv[int(num1[i]*2-2)]-miu*uv[int(num1[i]*2-1)]+lamda*uv[int(num2[i]*2-2)]+miu*uv[int(num2[i]*2-1)])/length[i]
            #print('%d%d号杆的力是'%(num1[i],num2[i]),bar_force,end='\n')
            f.write('%d%d号杆的力是'%(num1[i],num2[i]))
            f.write(str(bar_force))
            f.write('\n')
        f.write('\n')
        f.close()

def for_num(str_):
    try:
        return float(eval(str_))
    except:
        return sympy.Symbol(str_)

if __name__=='__main__':
    main()

qq_45885422

关注

2
点赞
踩
37

收藏

觉得还不错? 一键收藏
6
评论
基于Python的结构力学位移法编程求解

位移法的基本思路使以结点位移(广义位移)作为基本未知量，写出由未知位移表示的应变，由物理方程写出仍由未知位移表示的应力表达式，最后由平衡条件解出所有的未知位移。本代码可用于平面桁架结构的计算。
复制链接

扫一扫