########################################################################## # program: mdof_arbit_force_rk4.py # author: Tom Irvine # version: 1.4 # date: May 1, 2012 # description: Calculate the response of an MDOF system to an applied # force using the Runge-Kutta fourth order method. # # The input file must have two columns: time(sec) & force ########################################################################## from tompy import read_two_columns from tompy import GetInteger2 from tompy import read_array,enter_int,enter_float from ode_plots import ode_plots,ode_plots_sdof from generalized_eigen import generalized_eigen from numpy import zeros,dot,linspace,interp,pi from numpy import atleast_1d import matplotlib.pyplot as plt ######################################################################## print(" ") print(" Select units: ") print(" 1=English 2=metric ") iu=GetInteger2() print(" ") if(iu==1): print(" units: mass (lbm)") print(" stiffness (lbf/in)") print(" velocity (in/sec)") print(" displacement (in)") else: print(" units: mass (kg)") print(" stiffness (N/m)") print(" velocity (m/sec)") print(" displacement (m)") print(" ") M=read_array("mass matrix") K=read_array("stiffness matrix") nlength=len(atleast_1d(M)) if(iu==1): M/=386 print(" ") ndof,fn,MS =generalized_eigen(K,M,2) print " " print " Select damping input method " print " 1=uniform damping ratio" print " 2=damping ratio vector " idm=GetInteger2() if(idm==1): damp=zeros(ndof,float) print " " print " Enter uniform damping ratio " udamp = enter_float() for i in range(0,ndof): damp[i]=udamp else: damp=read_array(" Enter damping ratio vector") num=ndof max_fn=max(fn) min_fn=min(fn) if(min_fn>0): period_max=1/min_fn else: period_min=1/max_fn period_min=1/max_fn dtr=period_min/20 ######################################################################## v0=read_array("initial velocity vector") d0=read_array("initial displacement vector") ######################################################################## if(iu==1): print "Each force file must have two columns: time(sec) & force(lbf)" else: print "Each force file must two columns: time(sec) & force(N)" print " " print " Enter the number of force files " nff=enter_int() MAX=100000 tt=zeros((MAX,nff),float) ff=zeros((MAX,nff),float) force_dof=zeros(ndof,int) for i in range (0,ndof): force_dof[i]=-999 print " " print " Note: the first dof is 1 " for i in range (0,nff): ii=i+1 print " " print " Enter force file %d " % ii a,b,num=read_two_columns() L=len(a) if(L>MAX): L=MAX tt[0:L,i]=a[0:L] ff[0:L,i]=b[0:L] print " " print " Enter the number of dofs at which this force is applied" nfa=enter_int() for j in range (0,nfa): print " " if(j==0 and nfa==1): print " Enter the dof number for this force " if(j==0 and nfa>1): print " Enter the first dof number for this force " if(j>0 and nfa>1): print " Enter the next dof number for this force " nn=enter_int() nn=nn-1 force_dof[nn]=i #******************************************************************************* print " " print " Enter the duration (sec) " dur=enter_float() srr=20*max_fn print " " print " Enter the sample rate (recommend > %8.4g) " % srr sr=enter_float() dt=1/sr nt=int(dur/dt) t = linspace(0., dur, nt) #******************************************************************************* # # interpolate force # print " begin interpolation" FFI=zeros((nt,nff),float) for i in range (0,nff): tstart=tt[0,i] tt[:,i]=tt[:,i]-tstart # print max(t),max(tt[:,i]) last=MAX for j in range (1,MAX): if(tt[j,i]<=tt[(j-1,i)]): last=j break # print last tint=tt[0:last,i] fint=ff[0:last,i] FFI[:,i] = interp(t,tint,fint) print " end interpolation" #******************************************************************************* # nt_m1=int(nt-1) omegan=zeros(ndof,float) cc=zeros(ndof,float) kk=zeros(ndof,float) if(nlength>1): x=zeros((nt,ndof),float) y=zeros((nt,ndof),float) a=zeros((nt,ndof),float) nx=zeros((nt,ndof),float) ny=zeros((nt,ndof),float) na=zeros((nt,ndof),float) else: x=zeros(nt,float) y=zeros(nt,float) a=zeros(nt,float) nx=zeros(nt,float) ny=zeros(nt,float) na=zeros(nt,float) # transform initial conditions if(nlength==1): nd0=MS*(M*d0) nv0=MS*(M*v0) else: nd0=dot(MS.T,dot(M,d0)) nv0=dot(MS.T,dot(M,v0)) for j in range (0,ndof): if(nlength>1): nx[0,j]=nd0[j] ny[0,j]=nv0[j] omegan[j]=2*pi*fn[j].real cc[j]=2*damp[j]*omegan[j] kk[j]=(omegan[j])**2 else: omegan=2*pi*fn cc=2*damp*omegan kk=omegan**2 nx[0]=nd0 ny[0]=nv0 FA=zeros(ndof,float) FB=zeros(ndof,float) FF1=zeros(ndof,float) FF2=zeros(ndof,float) FF3=zeros(ndof,float) FF4=zeros(ndof,float) nFF1=zeros(ndof,float) nFF2=zeros(ndof,float) nFF3=zeros(ndof,float) nFF4=zeros(ndof,float) for i in range (0,nt_m1): ta=t[i] tb=t[i+1] h=tb-ta hd2=h/2 for j in range (0,ndof): j_index=force_dof[j] if(j_index!=-999): FA[j]=FFI[i,j_index] FB[j]=FFI[i+1,j_index] FF1= FA FF2= (FB + FA)/2. FF3= FF2 FF4= FB nFF1=dot(MS.T,FF1) nFF2=dot(MS.T,FF2) nFF3=dot(MS.T,FF3) nFF4=dot(MS.T,FF4) if(nlength>1): for j in range (0,ndof): c=cc[j] k=kk[j] X1= nx[i,j] Y1= ny[i,j] A1= nFF1[j]-c*Y1 -k*X1 X2= X1+Y1*hd2 Y2= Y1+A1*hd2 A2= nFF2[j]-c*Y2 -k*X2 X3= X1+Y2*hd2 Y3= Y1+A2*hd2 A3= nFF3[j]-c*Y3 -k*X3 X4= X1+Y3*h Y4= Y1+A3*h A4= nFF4[j]-c*Y4 -k*X4 nx[i+1,j]= X1+(h/6.)*( Y1 + 2*Y2 + 2*Y3 + Y4 ) # displacement ny[i+1,j]= Y1+(h/6.)*( A1 + 2*A2 + 2*A3 + A4 ) # velocity na[i+1,j]= nFF4[j] -c*ny[i+1,j] -k*nx[i+1,j] else: c=cc k=kk X1= nx[i] Y1= ny[i] A1= nFF1-c*Y1 -k*X1 X2= X1+Y1*hd2 Y2= Y1+A1*hd2 A2= nFF2-c*Y2 -k*X2 X3= X1+Y2*hd2 Y3= Y1+A2*hd2 A3= nFF3-c*Y3 -k*X3 X4= X1+Y3*h Y4= Y1+A3*h A4= nFF4-c*Y4 -k*X4 nx[i+1]= X1+(h/6.)*( Y1 + 2*Y2 + 2*Y3 + Y4 ) # displacement ny[i+1]= Y1+(h/6.)*( A1 + 2*A2 + 2*A3 + A4 ) # velocity na[i+1]= nFF4 -c*ny[i+1] -k*nx[i+1] for i in range (0,nt): if(nlength>1): X=nx[i,:] Y=ny[i,:] AA=na[i,:] x[i,:]=dot(MS,X) y[i,:]=dot(MS,Y) a[i,:]=dot(MS,AA) else: X=nx[i] Y=ny[i] AA=na[i] x[i]=(MS*X) y[i]=(MS*Y) a[i]=(MS*AA) #******************************************************************************* if(iu==1): a=a/386 else: a=a/9.81 fignum=1 if(nlength>1): ode_plots(fignum,nt,ndof,t,x,y,a,iu) else: ode_plots_sdof(fignum,nt,ndof,t,x,y,a,iu) plt.show()