夏亚-雷鸣
范文一：用叠加法求挠度和转角
当材料在线弹性范围内工作时，梁的挠度、转角均与载荷成线性关系．而且弯曲变形是很小的．因此，当梁上同时作用几种载荷时，任一载荷引起的变形，不会受到其他载荷的影响，即每种载荷对弯曲变形的影响是各自独立的。所以，几种载荷同时作用下梁的挠度和转角，等于各种载荷单独作用下挠度和转角的代数和，这就是求解弯曲变形的叠加法．当只需确定某些指定截面的挠度和转角时，应用叠加法是比较方便的．下面举例说明．
例 7-3 图 7-8 所示简支梁，承受均布载荷 q 和集中力偶 M 0 作用，已知 M 0 =ql2 。试求跨度中点的挠度 f c 和 A 截面的转角 θA 。
解：利用叠加法求解时，首先将 q , M0 同时作用下的简支梁 ( 图 7 -8a ) ，分解为 q 作用下的简支梁 ( 图 7-8b) 和 M 0 作用下的简支梁 ( 图 7 -8c ) ，然后，由表 7.1 查取结果叠加。
从表的第 9 栏查得均布载荷 q 作用下的中点挠度和 A 端面转角分别为
由表 7.1 第 5 栏查得集中力偶 M 0 作用下的中点挠度和 A 端面转角分别为
叠加以上结果，求得 q , M0 同时作用下的中点挠度和 A 截面转角为
f c 为负值，表示挠度向下．θA 为负值，表示 A 截面顺时针转动．
例 7-4 简支梁如图 7 — 10a 所示，在 2a 的长度上对称地作用有均布载荷 q. 试求梁中点挠度和梁端面的转角．
解：利用叠加法求解。由于简支梁上的载荷对跨度中点 C 对称，故 C 截面的转角应为零．因而从 C 截面取出梁的一半，可将其简化为悬臂梁，如图 7 — 10b 所示。梁上作用有均布载荷 q 和支座 B 的反力 R B = qa.这样，悬臂梁上 B 端面的挠度在数值上等于原梁中点 C 的挠度，但符号相反， B 端面的转角即为原梁 B 端面的转角．经这样处理后，应用叠加原理求解比较方便．
由表 7 · 1 的第 2 栏查得，当集中力 R B (=qa) 作用时 ( 图 7 — 10c ) ， B 端面的转角和挠度分别为
由表 7 · 1 的第 4 栏查得，当均布载荷 q 作用时 ( 图 7 — 10d) ， E 截面的转角和挠度分别为
由于 EB 梁段上无载荷作用，所以 q 引起 B 点的转角和挠度分别为
=
=
叠加上述结果，可得 B 端面的转角和挠度分别为
于是，原梁 ( 图 7 — 10a ) 中点 C 的挠度 f c 为
例 7-6 某一变截面外伸梁如图 7 — 11a 所示． AB 、 BC 段的抗弯刚度分别为 EI 1 和 EI 2 ，在 C 端面处受集中力 P 作用，求 C 端面的挠度和转角．
解：由于外伸梁是变截面的，故不能直接应用表 7 ． 1 中的结果．为此，必须将外伸梁分为 AB 、 BC 两段来研究．首先假设梁的外伸段 BC 是刚性的，研究由于简支梁 AB 的变形所引起的 C 截面的挠度和转角．然后，再考虑由于外伸段 BC 的变形所引起的 C 截面的挠度和转角．最后将其两部分叠加，得 C 截面的实际变形．
由于假设 BC 段为刚性，故可将 P 力向简支梁 AB 的 B 端简化，得 P 和 Pa ．P 力可由 B 支座的反力平衡，不会引起简支梁的弯曲变形。集中力偶 Pa 引起 B 截面的转角 ( 图 7 — 11 b) 由表 6 ． 1 查得
它引起 C 截面的转角和挠度分别为
在考虑 BC 段的变形时，可将其看作悬臂梁 ( 图 7 — 11c ) ，由表 6 · 1 查得，在 P 力作用下 C 截面的转角和挠角分别为
将图 7 — 11b 、 c 中的变形叠加后，求得 C 端面实际的转角和挠度分别为
例 7-7 在悬臂梁 AB 上作用线性分布载荷，如图 7-12 所示．试求自由端 B 点的挠度．
解：本例同样可以应用叠加法求解．将图中 dx 微段上载荷 qdx 看作集中力，查表 7 · 1 的第 3 栏求得微段载荷 qdx 作用下自由端 B 截面的挠度为
根据题意，线性分布载荷的表达式为
(1)
(2)
按照叠加原理，自由端 B 点的挠度应为 df B 的积分．将 (2) 式代入 (1) 式，积分得
f B 为负号，表示方向向下．
范文二：用积分求悬臂梁自由端的转角和挠度
8-1、用积分求悬臂梁自由端的转角和挠度。
解:弯矩方程: q A B 21 M(x),,q(l,x)EI 2L
微分方程:
''21 EIy,q(l,x)z2
1'223EIy,0.5qlx,0.5qlx，qx，Cz6积分求解: 112234EIy,qlx,qlx，qx，Cx，D0.25z624
''由边界条件: 得:C=0，D=0 x,0,,,y,0;y,0AAA
34qlql,y则: ,,,BB6EI8EIzz
8-7 试用叠加法求自由端截面的转角和挠度。
解:分解荷载为单独作用在结构上。则:
F=qL/6 q 3322qlqlFlFl,,,，,,,，，, cbqcFA B 24EI8EI6EI144EIEI C
L L/2
q yCq y,y，y,l/2，l/2，y,,1ccqcFbqbFcFθBq
A 323B EI C qllFllql,,,，,，,024EI26EI224EIL L/2
8-18、一变截面梁，求C截面的转角和挠度。
F 解:用叠加法计算。 2EI EI A
1、分解荷载如图 C B
a a
222F FaFaaFa5Fa,,,,，,,，，, CBCFa 14EI2EI2EI4EI2EI A C
y32223C1 B FaFaaFaFaFayyy()a,，,，，，，a a 12ccc32EI22EI22EI2EI3EI,,,F 33FaEI ,2EIYC2 C B a a
8-23、已知横梁的抗弯刚度EI，
C 1、 位移条件:y=Δ BBC
43EA a qlNlNay2、 ,,,,,q BBCA B EIEIEA83
EI L 43Al3、 Nq3，8(Al3aI)
C
B
N q A B
范文三：103用积分求悬臂梁自由端的转角和挠度
8-1、用积分求悬臂梁自由端的转角和挠度。
解：弯矩方程： q A B 21 M(x),,q(l,x)EI 2L
微分方程：
''21 EIy,q(l,x)z2
1'223EIy,qlx,qlx，qx，C0.50.5z6积分求解：
112234EIy,qlx,qlx，qx，Cx，D0.25z624
''由边界条件： 得：C=0，D=0 x,0,,,y,0;y,0AAA
34qlql则： ,y,,,BB6EI8EIzz
8-7 试用叠加法求自由端截面的转角和挠度。
F=qL/6 解：分解荷载为单独作用在结构上。则：q
A B 3322EI C qlqlFlFl,,,，,,,，，, cbqcFL L/2 24EI8EI6EI144EI
q
A B EI C yCq y,y，y,l/2，l/2，y,,θ1ccqcFbqbFcFBq L L/2 323qllFllql,,,，,，,024EI26EI224EI
8-18、一变截面梁，求C截面的转角和挠度。
F 解：用叠加法计算。
2EI EI A 1、分解荷载如图 C B
a a
222F FaFaaFa5Fa,,,,，,,，，, Fa CBC12EI 4EI2EI2EI4EIA C
B yC1 32223FaFaaFaFaFaa a yyy()a,，,，，，，12ccc32EI22EI22EI2EI3EI,,,
F 33Fa,EI 2EIYC2 C B a a
8-23、已知横梁的抗弯刚度EI，
C 1、 位移条件：y=Δ BBC
43EA a qlNlNa,,,,,y2、 q BBCA B 8EI3EIEA
EI L 43Al3、 Nq38(Al，3aI)
C
B
N q A B
范文四：8-1、用积分求悬臂梁自由端的转角和挠度。
8-1、用积分求悬臂梁自由端的转角和挠度。 解:弯矩方程: q A B 21M(x),,q(l,x) 2EI L
微分方程:
''21EIy,q(l,x) z2
1'2230.50.5EIy,qlx,qlx，qx，Cz6积分求解: 112234EIy,qlx,qlx，qx，Cx，D0.25z624
''x,0,,,y,0;y,0由边界条件: 得:C=0，D=0 AAA
34qlql则:,y ,,,BB6EI8EIzz
8-7 试用叠加法求自由端截面的转角和挠度。
解:分解荷载为单独作用在结构上。则:
F=qL/6 q 3322qlqlFlFl,,,，,,,，，, cbqcFA 24EI8EI6EI144EIB EI C
L L/2
q yCq y,y，y,l/2，l/2，y,,θ1ccqcFbqbFcFBq
A 323B EI C qllFllql,,,，,，,024EI26EI224EIL L/2
8-18、一变截面梁，求C截面的转角和挠度。
F 解:用叠加法计算。 2EI EI A
1、分解荷载如图 C B
a a
222F FaFaaFa5Fa,,,,，,,，，, CBC1Fa 4EI2EI2EI4EI2EI A C
yC1 32223B FaFaaFaFaFayyy()a,，,，，，，a a 12ccc32EI22EI22EI2EI3EI,,,F 33FaEI ,2EIYC2 C B a a
8-23、已知横梁的抗弯刚度EI，
C 1、 位移条件:y=Δ BBC
43EA a qlNlNay,,,,,2、 q BBCA B EIEIEA83
EI 4L 3AlNq3、 38(Al3aI)，
C
B
N q A B
范文五：用模态叠加法求固有频率
用模态叠加法求固有频率
一、 模态分析法(振型叠加法)原理
对于n个自由度系统，其在广义坐标系下的运动微分方程为
，， (1-1) MxkxFt，,(),,,,,,,,,,
设在t=0时，有初始条件:
和 ，，xx(0),xx(0),,,,,,,,,00
通过求解特征值问题，可得系统的固有频率和振型向量
,,(1,2,,)uin,？,,nii
1 ,,,？uin(1,2,,),,,,iiMi
以正则振型矩阵,作为变换矩阵，令 ,,
xz,, (a) ,,,,,,
T，，代入方程(1-1)，并前乘以正则振型矩阵的转置，得 ,，，
TTT，，,,,,,MzkzFt，,() (b) ,,,,,,,,,,,,,,,,,,,,
T,,,MI? ,,,,,,,,
2，，,n1,,2,,T, ,,,k,,,,,,,,,,n2,,？,,2,,,nn，，
TPtFt()(),,令 ---- 是正则坐标系下的激励。 ,,,,,,
则方程(b)为
，，zzPt，,,() (c) ,,,,,,,,
展开后，得
2,，，zzPt，,,()1111n,2,，，zzPt，,(), (1-2) 2222n,？？,2,，，zzPt，,(),nnnnn,
1
T式中 ，为对应第i个正则坐标的激励。 PtFtin()()(1,2,,),,,？,,,,ii
对于方程(1-2)是一组n个独立的方程，每个方程和单自由度系统的强迫振动相同，因此可按单自由度系统中的方法独立地求解每个方程。则由杜哈美积分得方程(1-2)的通解
，zi0,，,,ztztt()cossiniinini0,ni 1t，，,，,,？Ptdin()sin()1,2,,,,,,0nii,ni
式中和 是第i个正则坐标的初始位移和初始速度。 ，zzi0i0
? xz,,,,,,,,
? (d) xz,,,,,,,,00
和 (e) ，xz，,,,,,,,,00
T,M 前乘以式(d)两端，得 用,,,,
TTMxMz ,,,,,,,,,,,,,,,,,,00
T? zMx,,,,,,,,,,00
T，，同理，有 zMx,,,,,,,,,,00
写成分量形式
T zMxin,,,,(1,2,,)？,,,,,,i0i0
T，，zMxin,,,,(1,2,,)？ ,,,,,,i0i0
zx最后，由方程(a)，将正则坐标的解变换到原广义坐标，就得到方程(1-1)的,,,,
解。
(2)在许多工程问题中系统的自由度很多，要想求出系统的所有固有频率和振型向量，计算成本很大，有时甚至是不可能的。由于激励的高频成分很微弱，或者由于系统的高频振动没有激发出来，总之系统的响应中只有较低的几阶振型分量。因此，使用振型叠加法可以使计算大大地简化。
例如，若系统为n自由度，且只需考虑前p(p<>
2
p阶固有频率和振型向量:
,,,(1,2,,)ip,？,,nii
，，则振型矩阵 是一个n×p的矩阵。作如下的线性变换 ,,,,,？,,,,,,,,12p,,，，
xz,,,,,,,,np11，，np，
T,代入方程(1-1)，并前乘以正则振型矩阵的转置，得 ,,
TTT，，,,,,,MzkzFt，,() (b) ,,,,,,,,,,,,,,,,,,,,
T,,,MI ,,,,,,,,，，，，pnnnnppp
2，，,n1,,2,,T, n2k,,,,,,,,,,,,,,,pnnnnppp，，，，？,,2,,,pp，，pp，
TPtFt()(),, ,,,,,,pn，，11pn，
，，? zzPt，,,(),,,,,,,,ppp，，，111pp，
展开后，得
2,，，zzPt，,,()1111n,2,，，zzPt，,(),2222n ,？？,2,，，zzPt，,(),ppppp,
则原广义坐标下的响应为
xzzzz,,，，，,,,,？ ,,,,,,,,,,,,12pnpp1112，，np，
这种处理方法称为振型截断法。
二、振型叠加法计算步骤 1. 建立广义坐标系下的运动微分方程
，，MxkxFt，,() ,,,,,,,,,,
，给出广义坐标系下初始条件:xx和 ,,,,002. 计算固有频率和振型向量
3
,,(1,2,,)uin,？,,nii
T计算主质量: MuMuina,,(1,2,,)？,,,,,,iii
1得正则振型向量: ,,,？uin(1,2,,),,,,iiMi
3. 初始条件变换
T zMxin,,,(1,2,,)？,,,,,,i0i0
T， ，zMxin,,,(1,2,,)？,,,,,,i0i0
4. 激励变换
T PtFtin()()(1,2,,),,,？,,,,ii
5. 计算正则坐标系下的响应
正则坐标系下运动微分方程
2，，，,,,？ zzPtin()(1,2,,)iniii
利用杜哈美积分写出正则坐标系下的响应
，zi0,，,,ztztt()cossiniinini0,ni 1t,，,,？Ptdin()sin()(1,2,,),,,,0nii,ni
6. 写出原广义坐标系下的响应
xzzzz,,，，，,,,,？,,,,,,,,,,,,12n12n
三、源程序
****************************** * Program to generate the response of an N-DOF structure
* to any arbitrary applied loading using the * MODE SUPERPOSITION METHOD ******************************
* Program Author : Sanjoy Chakraborty * Auburn University ******************************
program mode_s
use msimsl
implicit real *8 (a-h,o-z)
real *8 mr,kr
4
character *12 flin,fout
dimension mr(10,1),cr(10,1),kr(10,1),phi(10,10),tt(20,10),
ft(20,10),fr(6000,10),q(6000,10),q1(6000,10),q2(6000,10),
qt(10,1),q1t(10,1),q2t(10,1),nfor(10),xmin(10,1),vmin(10,1)
data ndim/10/
*******************************
* Assign I/O Devices and read in system parameters *******************************
write(*,50)
50 format(///,' Program for Eigensolution and Mode Superposition'
,//,' Program Author : Sanjoy Chakraborty, Auburn University'
,//////,' Enter Input File Name : ')
read(*,'(a)')flin
open(1,file=flin,status='old')
write(*,'(//a)')' Enter Output File Name : '
read(*,'(a)')fout
open(5,file=fout,status='unknown')
call input1(mr,cr,kr,phi,tt,ft,nd,nfor,ndim,delt,ttot,xmin,vmin)
open(2,file='x.out',status='unknown')
open(3,file='x1.out',status='unknown')
open(4,file='x2.out',status='unknown')
*******************************
* Generate Modal Force Vectors at discrete time intervals * fr(i,j) contains the modal force vector Fr at time step 'i' * j=1,nd corresponds to the DOF at which fr is applied * {fr} = [Phi]'.{ft}
* where, {ft} = physical force vector at time step 'i' * [Phi] = modal transformation matrix for given system * [Phi]' = transpose of [Phi]
*******************************
call force1(tt,ft,delt,phi,nd,fr,nfor,ndim,ttot)
*******************************
* Generate MODAL time history for all DOF's
* q(i,j), q1(i,j), q2(i,j) contain the modal
* displacements, velocities and accelerations resp. * i --> corresponds to time step number
* j --> corresponds to a particular DOF
5
*******************************
* Direct Integration of the uncoupled equations of motion * are carried out using the Newmark-B method *******************************
do 10 i=1,nd
call newb(i,mr(i,1),cr(i,1),kr(i,1),fr,q,q1,q2,delt,ndim,ttot,
xmin(i,1),vmin(i,1)) 10 continue
*******************************
* Transform MODAL time histories back to PHYSICAL Coords. * {x(t)} = [Phi].{q(t)} ..etc.
* Print Output
*******************************
ic=1
do 20 t=0.,ttot,delt
do 30 i=1,nd
qt(i,1)=q(ic,i)
q1t(i,1)=q1(ic,i)
q2t(i,1)=q2(ic,i)
30 continue
call matmul(phi,qt,nd,mr,ndim)
write(2,43)t,(mr(nn,1),nn=1,nd) 43 format(f8.4,10e14.6)
do 40 i=1,nd
q(ic,i)=mr(i,1)
40 continue
call matmul(phi,q1t,nd,mr,ndim)
write(3,43)t,(mr(nn,1),nn=1,nd)
do 41 i=1,nd
q1(ic,i)=mr(i,1)
41 continue
call matmul(phi,q2t,nd,mr,ndim)
write(4,43)t,(mr(nn,1),nn=1,nd)
do 42 i=1,nd
q2(ic,i)=mr(i,1)
42 continue
ic=ic+1
6
20 continue
*************************
call res_spec(ic-1,q,q1,q2,delt,ndim,nd)
*************************
* Generate values of MAXIMUM RESPONSE
* @ each DOF
*************************
stop
end
*************************
subroutine input1(m,c,k,phi,t,f,nd,nfor,ndim,delt,ttot,xmin,
vmin)
*************************
* [m1] is the system mass matrix (input)
* [k1] is the system stiffness matrix (input)
* [m] contains the modal mass parameters (generated)
* [c] contains the modal damping parameters (generated) * [k] contains the modal stiffness parameters (generated) * alpha, beta = parameters to generate proportional damping matrix * [c] = alpha[m] + beta[k]
* [phi] contains the modal transformation matrix
* {t} contains the times at which the load vector
* {f} is specified
* nd = the number of degrees of freedom
* nsol = 1 for EIGENSOLUTION ONLY
* nsol = 2 for EIGENSOLUTION + MODE SUPERPOSITION ANALYSIS * delt = step size
* ttot = total time for which response will be evaluated
*************************
* This program uses the IMSL Math Subroutine DGVCSP to evaluate * the eigenvalues and eigenvectors of the dynamic system
*************************
implicit real *8 (a-h,o-z)
real *8 m,k,m1,k1
external dgvcsp
dimension m(ndim,1),c(ndim,1),k(ndim,1),phi(ndim,ndim),t(20,10)
dimension nfor(10),f(20,10),m1(10,10),c1(10,10)
dimension k1(10,10),omega(10)
dimension phit(10,10),eva(10),evec(10,10)
dimension temp1(10,10),temp2(10,10),temp3(10,10),xin(10,1),
7
vin(10,1),xmin(10,1),vmin(10,1)
read(1,*)nsol,nd
do 90005 i=1,nd
read(1,*)(m1(i,j),j=1,nd) 90005 continue
do 90006 i=1,nd
read(1,*)(k1(i,j),j=1,nd) 90006 continue
********************
do 825 i=1,nd
do 825 j=1,nd
temp1(i,j)=m1(i,j)
temp2(i,j)=k1(i,j) 825 continue
c Call the IMSL routine DGVCSP to evaluate the eigenvalues and c eigenvectors of the system.
call dgvcsp(nd,temp2,ndim,temp1,ndim,eva,evec,ndim)
do 903 i=1,nd
omega(i)=sqrt(eva(nd-i+1)) 903 continue
c write(*,*)(omega(i),i=1,nd)
c pause
c stop
do 904 i=1,nd
do 904 j=1,nd
phi(i,j)=evec(i,nd-j+1) 904 continue
******************************
* Normalize Eigenvectors wrt 1st row values ******************************
do 826 j=1,nd
temp=phi(1,j)
do 826 i=1,nd
8
phi(i,j)=phi(i,j)/temp 826 continue
******************************
* Print EIGEN OUTPUT
* Terminate if NSOL=1
******************************
write(5,905)
905 format(/,' E I G E N S O L U T I O N',//,
' Mode # Omega(rad/sec)',/,
' ----------------------------')
do 906 i=1,nd
write(5,907)i,omega(i) 907 format(i6,e22.8)
906 continue
write(5,912)
912 format(//,' Mode Shapes',/)
908 i=1,nd do
write(5,'(10E14.6)')(phi(i,j),j=1,nd) 908 continue
if(nsol.eq.1)stop
********************
********************
read(1,*)delt,ttot
read(1,*)alpha,beta
call trans(phi,phit,nd,ndim)
call matmul1(phit,m1,temp1,nd,ndim)
call matmul1(temp1,phi,temp2,nd,ndim)
do 909 i=1,nd
m(i,1)=temp2(i,i)
909 continue
call matmul1(phit,k1,temp1,nd,ndim)
call matmul1(temp1,phi,temp2,nd,ndim)
do 910 i=1,nd
k(i,1)=temp2(i,i)
910 continue
9
do 911 i=1,nd
zeta=0.5*(alpha/omega(i)+beta*omega(i))
c(i,1)=2.0*m(i,1)*omega(i)*zeta
911 continue
do 921 i=1,nd
do 921 j=1,nd
c1(i,j)=alpha*m1(i,j)+beta*k1(i,j)
921 continue
write(5,922)nd
922 format(//,' S Y S T E M P R O P E R T I E S',
//,' Degrees of Freedom =',i3,//,' Mass Matrix',/)
do 913 i=1,nd
write(5,914)(m1(i,j),j=1,nd) 914 format(10e15.6)
913 continue
write(5,915)
915 format(//,' Stiffness Matrix',/)
do 916 i=1,nd
write(5,914)(k1(i,j),j=1,nd) 916 continue
write(5,917)
917 format(//,' Damping Matrix',/)
do 918 i=1,nd
write(5,914)(c1(i,j),j=1,nd) 918 continue
write(5,923)(m(i,1),i=1,nd) 923 format(//,' Modal Masses',//,10e15.6)
write(5,924)(k(i,1),i=1,nd) 924 format(//,' Modal Stiffnesses',//,10e15.6)
write(5,925)(c(i,1),i=1,nd) 925 format(//,' Modal Damping',//,10e15.6)
read(1,*)(xin(i,1),i=1,nd)
read(1,*)(vin(i,1),i=1,nd)
10
do 851 i=1,nd
do 851 j=1,nd
temp1(i,j)=0.0
temp2(i,j)=0.0
temp3(i,j)=0.0
851 continue
do 852 i=1,nd
temp1(i,i)=1.0/m(i,1)
852 continue
call matmul1(temp1,phit,temp2,nd,ndim)
call matmul1(temp2,m1,temp3,nd,ndim)
call matmul(temp3,xin,nd,xmin,ndim)
call matmul(temp3,vin,nd,vmin,ndim)
*********************************
* Input of forcing Function
*********************************
do 90011 i=1,nd
nfor(i)=2
t(1,i)=0.0
t(2,i)=ttot+100.
f(1,i)=0.0
f(2,i)=0.0
90011 continue
90014 read(1,*,err=90012)i,nfor(i)
do 90013 j=1,nfor(i)
read(1,*)t(j,i),f(j,i)
90013 continue
goto 90014
90012 continue
return
end
**************************
subroutine force1(tt,f,dt,p,nd,fr,nfor,ndim,ttot)
**************************
implicit real *8 (a-h,o-z)
dimension tt(20,10),f(20,10),fr(6000,ndim),p(ndim,ndim),
11
f1(10,1),pt(10,10),nfor(ndim),ft(10,1)
do 90030 i=1,nd
do 90030 j=1,nd
pt(i,j)=p(j,i)
90030 continue
ic=1
do 90020 t=0.,ttot,dt
*******************************
* Interpolate to obtain value of forcing function Ft * at time = t
*******************************
do 90022 nf=1,nd
do 90024 n=1,nfor(nf)-1
if(t.ge.tt(n,nf).and.t.le.tt(n+1,nf))then
ft(nf,1)=f(n,nf)+(f(n+1,nf)-f(n,nf))*
(t-tt(n,nf))/(tt(n+1,nf)-tt(n,nf))
goto 90022
endif
90024 continue
90022 continue
*******************************
* Generate MODAL force vector at time = t *******************************
call matmul(pt,ft,nd,f1,ndim)
do 90026 ii=1,nd
fr(ic,ii)=f1(ii,1) 90026 continue
ic=ic+1
90020 continue
return
end
*******************************
subroutine newb(ii,m,c,k,f,x,x1,x2,del,ndim,ttot,xin,vin)
*******************************
* Direct Integration of the uncoupled equations of motion * using the NEWMARK-b method
******************************
12
implicit real *8 (a-h,o-z)
real *8 m,k,ks
dimension f(6000,ndim),x(6000,ndim),x1(6000,ndim),x2(6000,ndim)
******************************
* Assign Initial Conditions (@ t=0.0)
******************************
x(1,ii)=xin
x1(1,ii)=vin
x2(1,ii)=(f(1,ii)-c*x1(1,ii)-k*x(1,ii))/m
ks=k+2.0*c/del+4.0*m/del**2
******************************
* Obtain time history at discrete time intervals of del
******************************
ic=1
do 101 t=0.,ttot,del
dps=f(ic+1,ii)+m*(4.0*x(ic,ii)/del**2+4.0*x1(ic,ii)
/del+x2(ic,ii))+c*(2.0*x(ic,ii)/del+x1(ic,ii))
x(ic+1,ii)=dps/ks
x2(ic+1,ii)=(4.0/del**2)*(x(ic+1,ii)-x(ic,ii)-
del*x1(ic,ii))-x2(ic,ii)
x1(ic+1,ii)=(2.0/del)*(x(ic+1,ii)-x(ic,ii))-x1(ic,ii)
c write(5,105)ic,t,f(ic+1,ii),f(ic,ii),dps,dx
c105 format('**',i6,e17.7,/,5x,2f15.9,/,5x,2e20.10)
ic=ic+1
101 continue
return
end
*****************************
subroutine matmul(a,b,nd,c,ndim)
******************************
* Used to multiply two matrices
* [C] = [A].{B}
*****************************
implicit real *8 (a-h,o-z)
dimension a(ndim,ndim),b(ndim,1),c(ndim,1)
do 9030 i=1,nd
do 9030 j=1,1
13
c(i,j)=0.0
do 9030 ii=1,nd
c(i,j)=c(i,j)+a(i,ii)*b(ii,j)
9030 continue
return
end
*****************************
subroutine res_spec(n,x,x1,x2,del,ndim,nd)
*****************************
* Obtain the maximum value of the response * (displacement, velocity or acceleration) * at each DOF for the time span considered ****************************
implicit real *8 (a-h,o-z)
dimension x(6000,ndim),x1(6000,ndim),x2(6000,ndim)
dimension xl(10),x1l(10),x2l(10),t(10),t1(10),t2(10)
do 9070 i=1,nd
xl(i)=x(1,i)
x1l(i)=x1(1,i)
x2l(i)=x2(1,i)
t(i)=0.
t1(i)=0.
t2(i)=0.
9070 continue
do 9060 i=2,n
do 9060 j=1,nd
if(abs(x(i,j)).gt.abs(xl(j)))then
xl(j)=x(i,j)
t(j)=float(i-1)*del
endif
if(abs(x1(i,j)).gt.abs(x1l(j)))then
x1l(j)=x1(i,j)
t1(j)=float(i-1)*del
endif
if(abs(x2(i,j)).gt.abs(x2l(j)))then
14
x2l(j)=x2(i,j)
t2(j)=float(i-1)*del
endif
9060 continue
write(5,9068)
9068 format(///,' R E S P O N S E S P E C T R A',///,
' DOF Maxm. @ t Maxm. @ t ',
' Maxm. @ t',/,
' X Vel. ',
' Accln.',/,
' ----------------------------------------------------',
'----------------')
do 9065 i=1,nd
write(5,9066)i,xl(i),t(i),x1l(i),t1(i),x2l(i),t2(i)
9066 format(i5,e15.7,f7.3,e15.7,f7.3,e15.7,f7.3)
9065 continue
return
end
*************************************
subroutine matmul1(a,b,c,nd,ndim) *************************************
implicit real*8 (a-h,o-z)
dimension a(10,10),b(10,10),c(10,10)
do 80010 i=1,nd
do 80010 j=1,nd
c(i,j)=0.0
do 80010 k=1,nd
c(i,j)=c(i,j)+a(i,k)*b(k,j) 80010 continue
return
end
***************************************
subroutine trans(p,pt,nd,ndim) ***************************************
15
implicit real*8 (a-h,o-z)
dimension p(10,10),pt(10,10)
do 80001 i=1,nd
do 80001 j=1,nd
pt(i,j)=p(j,i) 80001 continue
return
end
***********************************************
16
转载请注明出处范文大全网 » 用叠加法求挠度和转角
