回复 15 # fawcgzmg 的帖子
解方程的方法和程序大多都是通用的，我倒是建议你给作者发个邮件，论文上的联系方式就是干这个使的。
回复 16 # Rainyboy 的帖子
日本论文，没办法
楼主太牛了，我从来没想过还能用python实现这个动力学问题。最近我也一直在学python，主要是用于ABAQUS的后处理，不知楼主对这方面有什么研究没有？
楼主您好！我学习python是为了将它与abaqus结合使用，前几天我编了一个从odb文件中提取时程曲线的代码：from odbAccess import *
from string import *
from abaqusConstants import *
odb=openOdb('F:/abaqus odb file/0.3gtianjin/tianjin3g.odb')
len1=len(odb.steps['Step-1'].frames)
len2=len(odb.steps['Step-2'].frames)
frameRepository1=odb.steps['Step-1'].frames
frameRepository2=odb.steps['Step-2'].frames
f1=open('d4.dat','w+')
f2=open('d5.dat','w+')
f3=open('d6.dat','w+')
for frame1 in range(len1):
dis=frameRepository1.fieldOutputs['U'].getSubset(NODAL).values
dataline1=dis.data
dataline2=dis.data
dataline3=dis.data
dataline4=dis.data
data1=dataline1-dataline2
data2=dataline2-dataline3
data3=dataline3-dataline4
f1.write(str(data1)+'\n')
f2.write(str(data2)+'\n')
f3.write(str(data3)+'\n')
for frame in range(len2):
dis=frameRepository2.fieldOutputs['U'].getSubset(NODAL).values
dataline1=dis.data
dataline2=dis.data
dataline3=dis.data
dataline4=dis.data
data1=dataline1-dataline2
data2=dataline2-dataline3
data3=dataline3-dataline4
f1.write(str(data1)+'\n')
f2.write(str(data2)+'\n')
f3.write(str(data3)+'\n')
f1.close()
f2.close()
f3.close()但这段代码运行的效率实在是有点低啊。麻烦您帮我看下可以如何提高它的效率，谢谢！
这个不错
路过，支持一下
加油。很不错。
Rainyboy 发表于 2011-9-21 15:34 static/image/common/back.gif
回复 15 # fawcgzmg 的帖子
解方程的方法和程序大多都是通用的，我倒是建议你给作者发个邮件，论文上的联系 ...
您好，我用matlab编程用wilson,newmark和解析解计算了单自由度和两自由度的位移响应，计算结果是完全吻合的，但是对于多自由度，我不知道怎么办了
%function = Wilson( K, M, C, f, d1, v1, dt, tend )
% 利用Wilson-theta (Newmark) 解析解法计算结构的动力响应
%K=sysK;% K ----- 调用刚度矩阵
%M=sysM; % M ----- 调用质量矩阵
% C=0.03*M+0.02*K;% C ----- 阻尼矩阵
K=;
M=;
C=0;
dt=0.1;% dt ----- 时间步长
tend=40;% tend --- 结束时间
w=600;
= size( K ) ;
%Wilson算法
theta = 1.37 ;
tao = theta*dt ;
alpha0 = 6/tao^2 ;
alpha1 = 3/tao ;
alpha2 = 2*alpha1 ;
alpha3 = tao/2 ;
alpha4 = alpha0/theta ;
alpha5 = -alpha2/theta ;
alpha6 = 1-3/theta ;
alpha7 = dt/2 ;
alpha8 = dt^2/6 ;
K1 = K + alpha0*M + alpha1*C ;
d = zeros( n, floor(tend/dt) + 1 ) ;
v = zeros( n, floor(tend/dt) + 1 ) ;
a = zeros( n, floor(tend/dt) + 1 ) ;
f = zeros( n, floor(tend/dt) + 1 ) ;
d(1,1) = 1.2 ; %初始位移
d(2,1)=0;
v(:,1) = 0 ;%初始速度
f(:,1) =0 ;%初始载荷
a(:,1) = inv(M)*(f(:,1)-K*d(:,1)-C*v(:,1)) ;%初始加速度
t=0:dt:tend;
for i=2:1:length(t)
%t(i) = (i-1)*dt ;
f(:,i)=0*100*sin(w*t(i));
ftheta = floor(theta) ;
%fq = f(i-1+ftheta-1)+ (theta-ftheta)*( f(i+ftheta-1) - f(i+ftheta-2) ) ;
fq=f(:,i-1)+ftheta*(f(:,i)-f(:,i-1));
f1 = fq + M*(alpha0*d(:,i-1)+alpha2*v(:,i-1)+2*a(:,i-1))+ C*(alpha1*d(:,i-1)+2*v(:,i-1)+alpha3*a(:,i-1)) ;
dq= inv(K1)*f1 ;
a(:,i) = alpha4*(dq-d(:,i-1)) + alpha5*v(:,i-1) + alpha6*a(:,i-1) ;
v(:,i) = v(:,i-1) + alpha7 * ( a(:,i) + a(:,i-1) ) ;
d(:,i) = d(:,i-1) + dt*v(:,i-1) + alpha8 * ( a(:,i)+2*a(:,i-1) ) ;
end
%解析解
for i=1:length(t)
x(i)=0.4*cos(0.5*t(i))+0.8*cos(1.581*0.5*t(i));
end
%Newmark算法
gama = 0.5 ;
beta = 0.25 ;
= size( K ) ;
Nalpha0 = 1/beta/dt^2 ;
Nalpha1 = gama/beta/dt ;
Nalpha2 = 1/beta/dt ;
Nalpha3 = 1/2/beta - 1 ;
Nalpha4 = gama/beta - 1 ;
Nalpha5 = dt/2*(gama/beta-2) ;
Nalpha6 = dt*(1-gama) ;
Nalpha7 = gama*dt ;
NK1 = K + Nalpha0*M + Nalpha1*C ;
Nd = zeros( n, floor(tend/dt) + 1 ) ;
Nv = zeros( n, floor(tend/dt) + 1 ) ;
Na = zeros( n, floor(tend/dt) + 1 ) ;
Nd(1,1) = 1.2 ; %初始位
Nd(2,1)=0;
f(:,1) =0 ;%初始载荷
Na(:,1) = inv(M)*(f(:,1)-K*Nd(:,1)-C*Nv(:,1)) ;%初始加速度
t=0:dt:tend;
for i=2:1:length(t)
f(:,i)=0*100*sin(w*t(i));
f2 = f(:,i) + M*(Nalpha0*Nd(:,i-1)+Nalpha2*Nv(:,i-1)+Nalpha3*Na(:,i-1))+ C*(Nalpha1*Nd(:,i-1)+Nalpha4*Nv(:,i-1)+Nalpha5*Na(:,i-1)) ;
Nd(:,i) = inv(NK1)*f2 ;
Na(:,i) = Nalpha0*(Nd(:,i)-Nd(:,i-1)) - Nalpha2*Nv(:,i-1) - Nalpha3*Na(:,i-1) ;
Nv(:,i) = Nv(:,i-1) + Nalpha6*Na(:,i-1) + Nalpha7*Na(:,i) ;
end
plot(t,d(1,:),'-b^',t,Nd(1,:),'g+',t,x,'r')
xlabel('t');
ylabel('u(t)');
title('位移与时间的关系');
ME! 发表于 2012-12-19 15:47 static/image/common/back.gif
您好，我用matlab编程用wilson,newmark和解析解计算了单自由度和两自由度的位移响应，计算结果是完全吻合 ...
对于多自由度就不要固定矩阵的维数，编程时用一个变量来处理矩阵的维数呗
不是很懂您的意思，能说的详细点吗？怎么实现啊
本帖最后由 Rainyboy 于 2012-12-22 14:45 编辑
ME! 发表于 2012-12-22 09:47 static/image/common/back.gif
不是很懂您的意思，能说的详细点吗？怎么实现啊
以单步欧拉方法为例，在求解当前时间步的响应的程序中可以用矩阵运算，这样就不用关系模型的维数是多少了。
#本文件包含：瞬态响应分析方法之一：欧拉法
import numpy as math;
from DynamicProcessBasic import *;
from MechanicMode import *;
class EularMethod(LinearDynamicMethod):
'''欧拉法'''
def solve(self,dt,timezone):
LinearDynamicMethod.solve(self,dt,timezone);
#以局部变量引用，增加代码的可读性
(TotalSteps,dis_alltime,vec_alltime,acc_alltime,t_alltime) = \
(self.TotalSteps,self.dis_alltime,self.vec_alltime,self.acc_alltime,self.t_alltime);
#引入初始条件
t_alltime =self.timestart;
dis_alltime =self.Mode.InitCondition[:,0].T;
vec_alltime =self.Mode.InitCondition[:,1].T;
acc_alltime =self.Mode.InitCondition[:,2].T;
#开始递推
for i in range(1,TotalSteps):
t_alltime = self.timestart + i*dt;
dis_alltime = dis_alltime + self.dt*vec_alltime;
vec_alltime = vec_alltime + self.dt*acc_alltime;
acc_alltime = (self.Mode.M.I*(-1.0*self.Mode.K*(dis_alltime.T)- \
1.0*self.Mode.C*(vec_alltime.T)+(self.Mode.GetForceVectorAt(t_alltime)).T)).T;
我是这样做的，我是想说为什么我代入二自由度的质量矩阵和刚度矩阵求解时，求解的结果是正确的，代入多自由度的问题时就不行，可能我编程用的matlab，而您用的是python，可能也有点出入！
ME! 发表于 2012-12-22 16:39 static/image/common/back.gif
我是这样做的，我是想说为什么我代入二自由度的质量矩阵和刚度矩阵求解时，求解的结果是正确的，代入多自由 ...
哈哈，这样啊，会错意了，我是用的PYTHON，MATLAB以前用过一些，后来就全换成PYTHON了，你的多自由度测试算例在哪本书上找到的？还是拿别的软件算的？
在徐荣桥的结构分析的有限元法与matlab程序设计
徐斌-Matlab有限元结构动力学分析与工程应用
弱弱地问一句：为什么不用MATLAB编？