如果想用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,codecsdef main():path=input('请输入文件名：')try:file=codecs.open(path,'r',encoding='utf-8')except:print('文件不存在，请检查，程序结束运行。\n')data=file.readlines()time=0line=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+1for 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+1for 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+1file.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 mdef get_location_in_list(x, target):step = -1items = list()for i in range(x.count(target)):y = x[step + 1:].index(target)step = step + y + 1items.append(step)return itemsdef add_add(point,dis,add,loca):temp=[0]*(point*2)t=0for i in range(point*2):if i not in loca:temp[i]=add[t,0]t=t+1else:temp[i]=0return tempdef 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()

基于Python的结构力学位移法编程求解相关推荐

基于Fortran的结构力学位移法编程求解
本代码可用于平面桁架结构的计算. 本文是根据基于Python的结构力学位移法编程求解的另一版本的更新,内容基本相同,背景介绍此处省略.仅展示最终代码的运行过程. 由于Fortran编程语言使用比较困难 ...
matlab 矩阵位移法编程结构力学,matlab 矩阵位移法编程结构力学
矩阵位移法编程大作业 (091210211) 一.编制原理本程序的原理是基于结构力学矩阵位移法原理,以结构结点位移作基本未知量,将要分析的结构拆成已知节点力-结点力位移关系的单跨梁集合,通过强令结构 ...
matlab 矩阵位移法编程结构力学,matlab 矩阵位移法编程结构力学.doc
matlab 矩阵位移法编程结构力学.doc 矩阵位移法编程大作业(091210211)一.编制原理本程序的原理是基于结构力学矩阵位移法原理,以结构结点位移作基本未知量,将要分析的结构拆成已知节点力 ...
基于 Python 的自然邻域法空间插值的实现与优化
接上期基于 Python 的自然邻域法空间插值的实现与思考. 上期说到,我们仅仅利用自然邻域法基础原理进行插值,会出现许多空值.异常值,且与ArcGIS相同分辨率.范围下的插值结果对比(对比图 ...
C语言练习题：三色球分组,编程计算三色球问题。若一个口袋中放有12个球，其中有3个红色的，3个白色的，6个黑色的，从中任取8个球，问共有多少种不同的颜色搭配？请用穷举法编程求解。
编程计算三色球问题.若一个口袋中放有12个球,其中有3个红色的,3个白色的,6个黑色的,从中任取8个球,问共有多少种不同的颜色搭配?请用穷举法编程求解. **输入格式:无 **输出格式:"i ...
pythonspark实践_基于Python的Spark Streaming Kafka编程实践
版权声明:本文为CSDN博主原创文章,未经博主允许不得转载. 说明 Spark Streaming的原理说明的文章很多,这里不做介绍.本文主要介绍使用Kafka作为数据源的编程模型,编码实践,以及一些 ...
基于 Python 的自然邻域法空间插值的实现与思考
自然邻域法是基于区域大小按比例对这些样本应用权重来进行插值 (Sibson 1981),该插值也称为 Sibson 或"区域占用 (area-stealing)"插值.其基本 ...
基于Python的GUI图形用户界面编程
[无限嚣张(菜菜)]:hello您好,我是菜菜,很高兴您能来访我的博客,我是一名爱好编程学习研究的菜菜,每天分享自己的学习,想法,博客来源与自己的学习项目以及编程中遇到问题的总结. 座右铭:尽人事,听 ...
矩阵位移法是用于求解杆系结构的计算机方法,结构力学的教学思路
结构力学的教学思路 2019-10-14 版权声明举报文章一.经典位移法(以下简称位移法)和矩阵位移法都是求解杆系结构的基本方法,是结构力学课程中两个十分重要的内容,两种方法都是结构力学课程中讲授 ...

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

目录