模拟景观悬臂桁架结构的解析模型

打开实时脚本

本例展示了如何在静域和频域找到一个普通悬臂桁架结构的关节位移的参数化解析表达式。静态情况的解析表达式是精确的。频率响应函数的表达式是实际频率响应的近似降阶形式。

这个例子使用了以下符号数学工具箱™功能:

Simscape™代码生成使用simscapeEquation
利用LU分解陆
MATLAB®函数生成使用matlabFunction

定义模型参数

本例的目标是分析地联系位移d为均匀横截面面积参数一个对于悬臂桁架结构中的某一特定杆。控制方程表示为 $米 \frac{d^{2} x}{d t^{2}} + K x ＝ F$ ．在这里,d这是由悬臂的右上关节处的载荷得出的结果。桁架在最左边的连接处固定在墙上。

桁架杆由密度均匀的线弹性材料制成。横截面面积一个蓝色突出显示的条(见图)的值可以变化。所有其他参数，包括其他杆的均匀横截面面积，都是固定的。利用小位移和线性位移假设，得到尖端的位移。

首先，确定桁架的固定参数。

指定桁架的长度和高度以及顶部水平桁架杆的数量。

L = 1.5;H = 0.25;N = 3;

指定桁架杆材料的密度和弹性模量。

Rhoval = 1e2;Eval = 1e7;

指定其他桁架的横截面面积。

fixed = 1;

接下来，定义特定桁架单元的局部刚度矩阵。

计算局部刚度矩阵k把它表示成一个符号函数。通过局部刚度矩阵，将桁架单元的力和位移联系起来。桁架单元的每个节点具有水平位移和垂直位移两个自由度。每个桁架元素的位移是一个列矩阵，对应于[x_node1,y_node1,x_node2,y_node2]。得到的刚度矩阵是一个4 × 4矩阵。

信谊一个Elt真正的G = [cos(t) sint 0 0;0 0cost sint];kk = A*E/l*[1 -1;-1 1];k = simplify(G'*kk*G);localstiff fn = symfun(k，[t,l,E])

local刚度fn (t, l, E) =
                  
                    $\begin{array}{l} （ \begin{array}{c} σ_{2} & σ_{1} & - σ_{2} & - σ_{1} \\ σ_{1} & σ_{3.} & - σ_{1} & - σ_{3.} \\ - σ_{2} & - σ_{1} & σ_{2} & σ_{1} \\ - σ_{1} & - σ_{3.} & σ_{1} & σ_{3.} \end{array} ） \\ 在哪里 \\ σ_{1} ＝ \frac{一个 E 罪 （ 2 t ）}{2 l} \\ σ_{2} ＝ \frac{一个 E {因为 （ t ）}^{2}}{l} \\ σ_{3.} ＝ \frac{一个 E {罪 （ t ）}^{2}}{l} \end{array}$

用类似的方法计算质量矩阵。

信谊ρmm = A*rho*l/6*[2 1;1 2];m = simplify(G'*mm*G);localMassFn = symfun(m，[t,l,rho]);

现在，定义桁架的杆为矩阵酒吧．杆件、起始节点和结束节点的指标如下图所示。

矩阵的行数酒吧是桁架的杆数。的列酒吧有四个元素，分别表示以下参数:

长度(l）
与横轴的夹角(t）
启动接头指标
端接头指标

bars = 0 (2*N+2*(N-1)，4);

定义上条和对角条。注意，我们感兴趣的柱是第(N+1)个柱或柱号4，它是左边第一个对角线柱。

为n = 1: n leelem = L/ n;Bars (n，:) = [leelem,0,n,n+1];lelem =√((L/N)²+H²);bars(N+ N，:) = [lelem,atan2(H,L/N)，N+1+ N,N+1];结束

定义下条和竖条。

为n = 1: n -1 leelem = L/ n;bars(2*N+ N，:) = [leelem,0,N+1+ N,N+1+ N+1];lelem = H;bars(2*N+N-1+ N，:) = [lelem,pi/2,N+1+ N+1,N+1];结束

将条形图组装成符号全局矩阵

桁架结构有7个节点。每个节点有两个自由度(水平位移和垂直位移)。这个系统的自由度总数是14。

numdfs = 2*2*(N+1)-2

numdfs = 14

系统质量矩阵米和系统刚度矩阵K是符号矩阵。将每个条形图的局部元素矩阵组合成相应的全局矩阵。

K = sym(零(numdfs, numdfs));M = sym(零(numdfs, numdfs));为J = 1:size(bars,1)为杆j构造刚度和质量矩阵。Lelem = bars(j,1);telm = bars(j,2);kelem = local刚度fn (telem, leelem,Eval);melem = localMassFn(telem,lelem,rhoval);N1 = bars(j,3);N2 = bars(j,4);对于不在参数化区域A内的条，设置刚度%和质量矩阵的数值值。注意，尽管值kelem和melem没有符号参数，他们仍然是%符号对象(符号数字)。如果j ~= N+1 kelem = subs(kelem,A, fixed);melem = subs(melem,A，固定);结束排列索引。Ix = [n1*2-1,n1*2,n2*2-1,n2*2];在系统全局矩阵中嵌入局部元素。K(ix,ix) = K(ix,ix) + kelem;M(ix,ix) = M(ix,ix) + melem;结束

消去附着在壁面上的自由度。

wallDofs = [1,2,2*(N+1)+1,2*(N+1)+2];%要消除的自由度K(wallDofs，:) = [];K(:，wallDofs) = [];M(walldfs，:) = [];M(:，wallDofs) = [];

F负载矢量是否带负号Y在最右上角的关节方向。

F = 0 (size(K,1)，1);F(2* n) = -1;

提取Y位移在最右上角的关节，创建一个选择向量。

c = 0 (1,size(K,1));c(2*N) = 1;

从静态情况下的精确符号解创建模拟场景方程

求位移的解析解并绘制d作为函数一个．在这里,K \ F表示所有关节和处的位移c选择特定的位移。首先，为了紧凑起见，使用16位数字表示数值的符号解决方案。

d_Static = simplify(c*(K\F));vpa (d_Static, 16)

ans =
                  
                    $- \frac{0.00000001118033988749895 (504.7737197247586 {一个}^{2} + 384.7213595499958 一个 + 22.3606797749979)}{一个 (1.28 一个 + 0.8944271909999158)}$

fplot (d_Static [AFixed / 10、10 * AFixed]);持有在；包含(的横截面,）;ylabel (的位移,d '）;标题(”）

图中包含一个轴对象。axis对象包含一个functionline类型的对象。

将结果转换为Simscape方程ssEq在Simscape块中使用。要显示生成的Simscape方程，请删除分号。

信谊dssEq = simscapeEquation(d,d_Static);

显示表达式的前120个字符ssEq．

strcat (ssEq (1:120),“……”）

ans = ' d = = (i - e (5.0) * (* 2.2 + 2 + e * cos (9.272952180016122 e 1) * 2.0 + 2 + i - e (5.0) * ^ 2 * 1.16 + 2 + i - e (5.0) * 1.0 + 1 + 9.2729521800 ^ 2 * cos(……”

比较数字和符号静态解决方案金宝搏官方网站

比较精确的解析解和数值解在一个范围内一个参数值。这些值形成一个序列AFixed来5 * AFixed增量为1。

paramvalues = fixed *(1:5)';d_NumericArray = 0 (size(AParamValues));为k=1: number (AParamValues) beginDim = (k-1)*size(k,1)+1;endDim = k*size(k,1);创建一个数值刚度矩阵并提取数值解。KNumeric = double(subs(K,A,AParamValues(K)));d_NumericArray(k) = c*(KNumeric\F);结束

计算符号值AParamValues．为此，调用潜艇函数，然后使用将结果转换为双精度数双．

d_SymArray = double(subs(d_Static,A,AParamValues));

在表格中可视化数据。

T = table(paramvalues,d_SymArray,d_NumericArray)

T =5×3表AParamValues d_SymArray d_NumericArray  ____________ ___________ ______________ 1 -4.5488 -4.5488 -4.6885 -4.6885 e-06 e-06 2 e-06 e-06 3 -4.5022 e-06 -4.5022 e-06 4 -4.4649 -4.4649 -4.4789 -4.4789 e-06 e-06 5 e-06 e-06

频率响应的近似参数符号解

频率响应的参数表示H (,)能否有效地快速检测参数的影响一个的每一个新值都不需要进行可能昂贵的数值计算一个．然而，为一个具有额外参数的大型系统找到一个精确的封闭形式的解通常是不可能的。因此，对于这样的系统，通常需要近似求解。本例将参数频响在零频率附近的近似如下:

使用可变精度算术，加快计算速度(vpa）.
求传递函数的Padé近似 $H （年代，一个）＝ c （ {年代}^{2} 米（一个） + K （一个））^{- 1} F$ 在频率附近S = 0用这个系列的前三个瞬间。其思想是，给定一个传递函数，您可以计算Padé近似矩为 $c （ - K （一个）^{- 1} 米（一个））^{j} K （一个）^{- 1} F$ ,在那里 $j \in ｛ 0 ， 2 ， 4 ， 6 ，．．．｝$ 对应于幂级数项 $｛ 1 ， {年代}^{2} ， {年代}^{4} ， {年代}^{6} ，．．．｝$ ．对于本例，使用computeHApprox(年代)是前三项的和。这种方法是降低传递函数阶数的一种非常基本的方法。
为了进一步加快计算速度，用的泰勒级数展开来近似每个矩项的内项一个周围AFixed．

使用vpa加快计算速度。

数字(32);K = vpa(K);M = vpa(M);

计算的LU分解K．

[LK,UK] = lu(K);

创建辅助变量和函数。

iK = @(x) UK\(LK\x);iK_M = @(x) -iK(M*x);momentPre = iK(F);

定义一个与前三个力矩对应的频率序列。在这里,年代是频率。

信谊年代sPowers = [1 s^2 s^4];

设定第一时刻，这是一样的d_Static在之前的计算中。

力矩= d_Static;

计算剩余的力矩。

为p=2: number (sPowers) momentPre = iK_M(momentPre);为pp=1: nummel(动量pre)动量pre (pp) = taylor(动量pre (pp)，A, fixed，“秩序”3);结束moment = c*momentPre;Moments = [Moments;moment];结束

结合力矩项来创建近似的解析频响Happrox．

HApprox = sPowers*moments;

显示前三个瞬间。因为表达式很大，所以只能部分显示它们。

arrayfun(@char,moments，“UniformOutput”、假);freqString = arrayfun(@char,sPowers.'，“UniformOutput”、假);表(freqString momentString,“VariableNames”, {“FreqPowers”，“时刻”})

ans =3×2表FreqPowers时刻  __________ _____________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________ {' 1 '}{”——(5 ^ (1/2)* (220 +200 * 5 * cos (8352332796509007/9007199254740992) + 116 * ^ (1/2) * ^ 2 + 10 * 5 ^ (1/2) + 40 * ^ 2 * cos (8352332796509007/9007199254740992) + 40 * ^ 2 + 20 * 5 ^ (1/2) * + 152 * 5 ^ (1/2) * ^ 2 * cos (8352332796509007/9007199254740992) + 36 * 5 ^ (1/2) * ^ 2 * cos (8352332796509007/4503599627370496))) / (200000000 * A * (A - A * cos (8352332796509007/4503599627370496) - 5 ^ (1/2) * cos(8352332796509007/9007199254740992) + 5 ^(1/2)))”}{' s ^ 2}{“0.000000000084654421099019119162064316556812 * (- 1)^ 2 - 0.000000000079001242991597407593795324876303 * +0．0000000004452470882909076005654298481594' } {'s^4'} {'0.000000000000012982169746567833536274593916742*A - 0.000000000000015007141035503798552656353179406*(A - 1)^2 - 0.000000000000045855285883001825767717087472451' }

将频率响应转换为包含数值的MATLAB函数一个而且年代．您可以在没有符号数学工具箱的情况下使用得到的MATLAB函数。

HApprox = matlabFunction“var”, [s]);

比较纯数字和符号派生的频率响应

改变频率0来1在logspace，然后将范围转换为弧度。

Freq = 2*pi*logspace(0,1);

计算频率响应的数值A =固定*摄动因子，即对于周围的小扰动一个．

摄动因子= 1 + 0.25;KFixed = double(subs(K,A, fixed *摄动因子));MFixed = double(subs(M,A, fixed *摄动因子));H_Numeric = 0 (size(freq));为k=1:数值(频率)sVal = 1i*频率(k);H_Numeric(k) = c*((sVal^2*MFixed + KFixed)\F);结束

计算频率响应在和频率范围内的符号值A =扰动因子*固定．

H_Symbolic = HApproxFun(固定*摄动因子，频率*1i);

画出结果。

图重对数(频率/(2 *π),abs (H_Symbolic),频率/(2 *π),abs (H_Numeric),“k *’）;包含(“频率”）;ylabel (的频率响应）;传奇(“象征”，“数字”）;

图中包含一个轴对象。axis对象包含2个line类型的对象。这些对象代表符号、数字。

分析曲线和数值曲线在所选区间内接近，但一般来说，在的邻域之外的频率曲线发散S = 0．同样，对于的值一个远离AFixed，泰勒近似引入更大的误差。