你好，请问一下桁架结构带未知参数的刚度矩阵如何在matlab中实现？

你提出的问题我之前刚好做过,使用有限元方法来进行桁架结构分析.

Matlab编程实现平面杆单元分析

首先,明确Matlab程序要实现的5个重要模块分别为：单元刚度矩阵的求解、单元组装、节点位移的求解、单元应力的求解、节点力的求解.下面给出这5个模块的实现.

1．\x09单元刚度矩阵求解

定义函数Bar2D2Node_Stiffness,该函数计算单元的刚度矩阵,输入d性模量E,横截面积A,两个节点坐标输出单元刚度矩阵k(4X4).具体代码如下：

function k=Bar2D2Node_Stiffness(E,A,x1,y1,x2,y2)

L=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1))

x=acos((x2-x1)/L)

C=cos(x)

S=sin(x)

k=E*A/L*[C*C C*S -C*C -C*SC*S S*S -C*S -S*S

-C*C -C*S C*C C*S-C*S -S*S C*S S*S]

2．\x09单元组装

定义函数Bar2D2Node_Assembly,该函数进行单元刚度矩阵的组装,输入单元刚度矩阵k,单元的节点编号i、j.输出整体刚度矩阵KK,具体代码如下：

function z = Bar2D2Node_Assembly(KK,k,i,j)

DOF(1)=2*i-1

DOF(2)=2*i

DOF(3)=2*j-1

DOF(4)=2*j

for n1=1:4

for n2=1:4

KK(DOF(n1),DOF(n2))= KK(DOF(n1),DOF(n2))+k(n1,n2)

end

end

z=KK

3．\x09节点位移的求解

定义函数Bar2D2Node_Disp(KK,num,p),该函数输入KK为总体刚度矩阵；num为活动自由度编号数组；p为活动自由度方向上的节点力；输出节点位移列阵.具体代码如下：

function u = Bar2D2Node_Disp(KK,num,p)

k=KK(num,num)

u=k\p

4．\x09单元应力的求解

定义函数函数Bar2D2Node_Stress(E,x1,y1,x2,y2,u),该函数计算单元的应力输入d性模量E,第一个节点坐标（x1,y1）,第二个节点坐标（x2,y2）单元节点位移矢量u,返回单元应力标量stress .具体代码如下：

L=sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1))

x=acos((x2-x1)/L)

C=cos(x)

S=sin(x)

stress=E/L*[-C -S C S]*u

5．\x09计算节点力

定义函数Bar2D2Node_Forces(KK,q),该函数用于计算节点力,KK为刚度矩阵,q为节点位移阵列

function P= Bar2D2Node_Forces(KK,q)

q=zeros(8,1)

q(num)=u

P=KK*q

至此,基于Matlab的杆单元有限元分析的程序设计已经完成,遇到实际问题时可以直接调用这些函数就可以解决问.

经典算例

如图所示的结构,各个杆的d性模量和横截面积都为 , .试基于MATLAB平台求解该结构的节点位移、单元应力以及支反力.

四杆桁架结构

对该问题进行有限元分析的过程如下

（1） 结构的离散化与编号

\x09对该结构进行自然离散,节点编号和单元编号如上图所示

（2）计算各单元的刚度矩阵(基于国际标准单位)

\x09\x09输入d性模量E、横截面积A,各点坐标.然后分别针对单元1,2,3和4,调用4次Bar2D2Node_Stiffness,就可以得到单元的刚度矩阵.

对应的主程序中代码：

E=2.95e11A=0.0001x1=0y1=0x2=0.4y2=0x3=0.4y3=0.3x4=0y4=0.3

k1=Bar2D2Node_Stiffness (E,A,x1,y1,x2,y2)

k2=Bar2D2Node_Stiffness (E,A,x2,y2,x3,y3)

k3=Bar2D2Node_Stiffness (E,A,x1,y1,x3,y3)

k4=Bar2D2Node_Stiffness (E,A,x4,y4,x3,y3)

（3） 建立整体刚度方程

由于该结构共有4个节点,因此,设置结构总的刚度矩阵为KK(8×8),先对KK清零,然后四次调用函数Bar2D2Node _Assembly进行刚度矩阵的组装.相关主程序代码为：

KK=zeros(8,8)

KK=Bar2D2Node_Assembly (KK,k1,1,2)

KK=Bar2D2Node_Assembly (KK,k2,2,3)

KK=Bar2D2Node_Assembly (KK,k3,1,3)

KK=Bar2D2Node_Assembly (KK,k4,4,3)

（4）边界条件的处理及刚度方程的求解

由图可以看出,被约束的自由度有：节点1的x,y方向自由度,节点2的y方向自由度,4节点的x、y方向两个自由度.则活动自由度编号为3,5,6.活动自由度对应的节点载荷F3=20000N,F5=0N,F6=25000N,采用高斯消去法进行求解,对应的代码为：

num=[3,5,6]%可活动的自由度编号

p=[200000-25000]

u=Bar2D2Node_Disp(KK,num,p)

（5）支反力的计算

在得到整个结构的节点位移后,由原整体刚度方程就可以计算出对应的支反力.这部分对应的主程序的代码如下：

q=zeros(8,1)

q(num)=u%节点位移阵列

P=Bar2D2Node_Forces(KK,q)

（6）单元应力的计算

先从整体位移列阵q中提取出单元的位移列阵,然后,调用计算单元应力的函数Bar2D2Node_ElementStress,就可以得到各个单元的应力分量.

u1=[q(1)q(2)q(3)q(4)]

stress1=Bar2D2Node_Stress(E,x1,y1,x2,y2,u1)

u2=[q(3)q(4)q(5)q(6)]

stress2=Bar2D2Node_Stress(E,x2,y2,x3,y3,u2)

u3=[q(1)q(2)q(5)q(6)]

stress3=Bar2D2Node_Stress(E,x1,y1,x3,y3,u3)

u4=[q(7)q(8)q(5)q(6)]

stress4=Bar2D2Node_Stress(E,x4,y4,x3,y3,u4)

（7）计算结果的整理

通过主程序的运行得计算结果.

主程序

%计算各单元的刚度矩阵(以国际标准单位)

E=2.95e11

A=0.0001

x1=0

y1=0

x2=0.4

y2=0

x3=0.4

y3=0.3

x4=0

y4=0.3

k1=Bar2D2Node_Stiffness (E,A,x1,y1,x2,y2)

k2=Bar2D2Node_Stiffness (E,A,x2,y2,x3,y3)

k3=Bar2D2Node_Stiffness (E,A,x1,y1,x3,y3)

k4=Bar2D2Node_Stiffness (E,A,x4,y4,x3,y3)

%建立整体刚度方程

%由于该结构共有4个节点,因此,结构总的刚度矩阵为KK(8×8),先对K清零,然后四次调用函数Bar2D2Node _Assembly进行刚度矩阵的组装.

KK=zeros(8,8)

KK=Bar2D2Node_Assembly (KK,k1,1,2)

KK=Bar2D2Node_Assembly (KK,k2,2,3)

KK=Bar2D2Node_Assembly (KK,k3,1,3)

KK=Bar2D2Node_Assembly (KK,k4,4,3)

%边界条件的处理及刚度方程求解

num=[3,5,6]%可活动的自由度编号

p=[200000-25000]

u=Bar2D2Node_Disp(KK,num,p)

%支反力的计算

q=zeros(8,1)

q(num)=u%节点位移阵列

P=Bar2D2Node_Forces(KK,q)

%各单元的应力计算

u1=[q(1)q(2)q(3)q(4)]

stress1=Bar2D2Node_Stress(E,x1,y1,x2,y2,u1)

u2=[q(3)q(4)q(5)q(6)]

stress2=Bar2D2Node_Stress(E,x2,y2,x3,y3,u2)

u3=[q(1)q(2)q(5)q(6)]

stress3=Bar2D2Node_Stress(E,x1,y1,x3,y3,u3)

u4=[q(7)q(8)q(5)q(6)]

stress4=Bar2D2Node_Stress(E,x4,y4,x3,y3,u4)

说的可能有些罗嗦,注意其中有5个function,和最后一个主程序,计算的时候直接运行主程序就可以了.希望能帮助到你.

两拉杆系统的刚度矩阵怎么求具体步骤如下：

1、弄清应力-应变关系以及刚度矩阵、柔度矩阵。

2、对于两拉杆系统材料，刚度矩阵可写为Qij。

3、Qij的计算公式可以通过Q11、Q12、Q22、Q23、Q66、∧计算可以得出。

4、复合材料工程常数之间的关系由麦克斯韦定理给出Vij。

5、刚度矩阵坐标转换公式为{CGij}。

6、应力张量转换矩阵即为刚度转换矩阵，与应变张量转换矩阵不同。

7、m，n，l为夹角余弦。（坐标轴1、2、3与坐标轴x、y、z夹角余弦分别为（li、mi、ni））。刚度矩阵就是节点载荷与节点位移之间的关系单元刚度矩阵xi单元e是在节点力作用下处于平衡。

机械振动刚度矩阵怎么求:

在 机械振动理论(3)-解析实模态分析 中，详细介绍了解析实模态分析方法，并以两自由度振动系统为例，给出了解析模态分析的一般步骤：

根据结构的几何形状、边界条件和材料属性构建质量矩阵、刚度矩阵、阻尼矩阵，进而得到系统的运动方程；

如果是比例阻尼，忽略阻尼项，得到系统特征多项式；

根据系统的特征多项式，求解系统的特征向量，并将各个特征向量代入得到各自对应的特征向量；

将特征向量以特征值的大小顺序按列排布，得到系统的特征向量矩阵（即，模态振型矩阵）；

通过坐标变换，使得质量矩阵、刚度矩阵、阻尼矩阵完全解耦，再运用单自由度的理论结果，结合转换到模态坐标下的初始条件，得到模态坐标下的位移解，再经过一次坐标变换，将模态坐标下的位移转变为物理坐标系下的解，也就是我们需要的时域振动解。

实验模态分析的一般步骤：

通过实验测量系统的频响函数；

从测得的频响函数来估计这些模态参数。

频响函数（Frequency Response Function，FRF）描述了频域中输出（通常采用加速度传感器测量）与输入（激励）之间的关系，在进行测量时，通常会用到两种激振（输入）方法，一种是力锤激励，另一种是激振器激励。但从理论角度讲，通过力锤激励法测试和激振器激励法测试所采集到的数据完全相同。

对称矩阵，这意味着 传递函数矩阵也是对称矩阵，满足Maxwell互易定理。

把这个时间域的矩阵变换到拉氏域（复变量为），并且假设初始位移和初始速度都为零（在描述受迫振动时，往往不太关心初始瞬态响应，在阻尼存在的情况下，瞬态响应会逐渐消失，所以直接令初始条件为零），则有：

或者式中 称为

动刚度矩阵。

由此也可以得到传递函数的定义：按照线性代数知识，一个矩阵的逆矩阵可以由其伴随矩阵计算出来：

根据 机械振动理论(2)-多自由度系统 中的一系列推导。

---