Simscape的悬臂桁架结构分析模型

打开直播脚本

该示例显示了如何找到用于在静态和频率域中的琐碎悬臂桁架结构的关节位移的参数化分析表达式。静态案例的分析表达是准确的。频率响应函数的表达式是实际频率响应的近似阶数。

此示例使用以下符号数学工具箱™功能：

Simscape代码生成使用simscapeequation.
鲁分解使用鲁
使用MATLAB生成函数matlabFunction

定义模型参数

该示例的目标是分析对位移有关D.到均匀的横截面区域参数一种针对悬臂桁架结构中某一特殊杆件。控制方程表示为 $m \frac{{D.}^{2} X}{D. {T.}^{2}} + K. X = F$ ．在这里,D.悬臂右上角的负载产生。桁架在最左边的关节处附着在墙上。

桁架杆由具有均匀密度的线性弹性材料制成。横截面区域一种在蓝色突出显示的棒（见图）可以变化。所有其他参数，包括其他条的均匀横截面区域是固定的。通过使用小型和线性位移假设获得尖端的位移。

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

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

l = 1.5;h = 0.25;n = 3;

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

Rhoval = 1E2;eval = 1E7;

指定其他桁架杆的横截面区域。

afixed = 1;

然后，定义特定桁架单元的局部刚度矩阵。

计算当地的刚度矩阵K.并将其表达为符号功能。桁架元素的力和位移通过局部刚度矩阵相关。桁架元素的每个节点具有两度自由度，水平和垂直位移。每个桁架元素的位移是一个列矩阵，其对应于[x_node1，y_node1，x_node2，y_node2]。得到的刚度基质是4×4的基质。

Syms.一种E.L.T.真实的g = [cos（t）sin（t）0 0;0 0 cos（t）sin（t）];kk = a * e / l * [1 -1; -1 1];k =简化（g'* kk * g）;localstiffnessfn = symfun（k，[t，l，e]）

localstiffnessfn（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. {COS. （ T. ）}^{2}}{L.} \\ σ_{3.} = \frac{一种 E. {罪 （ T. ）}^{2}}{L.} \end{array}$

以类似的方式计算质量矩阵。

Syms.rho.mm = a * rho * l / 6 * [2 1; 1 2];m =简化（g'* mm * g）;localmassfn = symfun（m，[t，l，rho]）;

现在，定义桁架的条作为矩阵酒吧．条形，起始关节和末端接头的索引如下图所示。

矩阵的行数酒吧是桁架的棒的数量。列的列酒吧有四个元素，它表示以下参数：

长度(L.）
相对于水平轴的角度（T.）
起始联合指数
终结指数

条=零（2 * n + 2 *（n-1），4）;

定义上条和对角线。注意，感兴趣的条是(N+1)第N条或第4条，这是从左边开始的第一个对角线条。

为了n = 1：n lelem = l / n;酒吧（n，:) = [lelem，0，n，n + 1];lelem = sqrt（（l / n）^ 2 + h ^ 2）;条（n + n，：）= [lem，atan2（h，l / n），n + 1 + n，n + 1];结尾

定义较低和垂直条。

为了lelem = L/ n;酒吧(2 * N + N:) = (lelem 0 N + 1 + N, N + 1 + N + 1];lelem = H;酒吧(2 * N + N - 1 + N:) = (lelemπ/ 2,N + 1 + N + 1, N + 1);结尾

将棒装配到象征性的全局矩阵中

桁架结构有七个节点。每个节点具有两度自由度（水平和垂直位移）。该系统中自由度的总数为14。

numDofs = 2 * 2 * 2 (N + 1)

numdofs = 14.

系统质量矩阵m和系统僵硬矩阵K.是符号矩阵。将本地元素矩阵从每个栏组装到相应的全局矩阵中。

K =符号(0 (numDofs numDofs));M =符号(0 (numDofs numDofs));为了j = 1:尺寸(酒吧、1)BAR J的％构建刚度和质量基质。Lelem =酒吧（J，1）;Telem =酒吧（J，2）;kelem = localstiffnessfn（Telem，Lelem，eval）;Melem = LocalMassfn（Telem，Lelem，Rhoval）;n1 =条（J，3）;n2 =条（J，4）;不在参数化区域A内的条形的百分比，设置刚度％和大规模矩阵到数值。请注意，虽然值％kelem和melem没有符号参数，它们仍然存在％符号对象（符号数字）。如果j〜= n + 1 kelem = subs（kelem，a，a，abixed）;Melem =潜艇（梅尔姆，A，Afix）;结尾%排列指数。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];％DOFs消除k（walldofs，:) = [];K（：，walldofs）= [];m（walldofs，:) = [];M（：，walldofs）= [];

F负载向量是否具有负载负载y直接接头的方向。

F = 0(大小(K, 1), 1);F (2 * N) = 1;

提取y右上角的位移，创建一个选择向量。

C =零（1，尺寸（k，1））;C（2 * n）= 1;

从精确的符号解决方案中创建SIMSCUPE方程，用于静态案例

找到并绘制位移的分析解决方案D.作为一个函数一种．在这里,k \ f.代表所有关节的位移和C选择特定的位移。首先，显示表示使用16位进行紧凑性的16位数字值的符号解决方案。

d_Static =简化(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]);持有上;Xlabel（'横截面，a'）;ylabel（'位移，d'）;标题(''）

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

Syms.D.ssEq = simscapeEquation (d, d_Static);

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

Strcat（SSEQ（1：120），'...'）

ans ='d ==（sqrt（5.0）*（a * 2.2e + 2 + a * cos（9.27295218001612222e-1）* 2.0e + 2 + sqrt（5.0）* a ^ 2 * 1.16e + 2 + sqrt（5.0）* 1.0e + 1 + a ^ 2 * cos（9.2729521800 ...'

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

比较一系列的确切分析解决方案和数字解决方案一种参数值。值表格来自序列ablix至5 *相反增量为1。

aparamvalues = afixed *（1：5）';d_numericarray = zeros（大小（aparamvalues））;为了k = 1：numel（aparamvalues）begindim =（k-1）*尺寸（k，1）+1;enddim = k *尺寸（k，1）;％创建数字刚度矩阵并提取数字解决方案。KNumeric =双(潜艇(K, AParamValues (K)));d_NumericArray (k) = c * (KNumeric \ F);结尾

计算符号值aparamvalues.．要做到这一点，就呼叫subs函数，然后将结果转换为双精度数双．

d_symarray = double（dim（d_static，a，aparamvalues））;

可视化表中的数据。

t =表（aparamvalues，d_symarray，d_numericArray）

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

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

频率响应的参数表示h（s，a）可以有效地快速检查参数的效果一种在不诉诸每个新价值的潜在昂贵的数字计算一种．但是，为大型系统找到具有附加参数的精确闭合液解决方案通常是不可能的。因此，您通常近似对这种系统的解决方案。此示例近似于零频率附近的参数频率响应，如下所示：

通过使用可变精度算术来加快计算（VPA.）。
找到转移函数的Padé近似 $H （ S. 那一种） = C （ {S.}^{2} m （一种） + K. （一种））^{- 1} F$ 周围的频率s = 0.用这个级数的前三个矩。这个想法是，给定一个传递函数，你可以计算Padé近似矩为 $C （ - K. （一种）^{- 1} m （一种））^{j} K. （一种）^{- 1} F$ ,在那里 $j \in {0. 那 2 那 4. 那 6. 那．．．}$ 对应于Power Series术语 ${1 那 {S.}^{2} 那 {S.}^{4.} 那 {S.}^{6.} 那．．．}$ ．对于此示例，计算Happox（s，a）作为前三个术语的总和。这种方法是减少传递函数顺序的非常基本的技术。
为了进一步加速计算，就泰勒序列扩展而言，近似每个时刻术语的内部期限一种大约ablix．

用VPA.加快计算。

数字（32）;k = VPA（k）;m = VPA（m）;

计算的LU分解K.．

[LK，UK] = Lu（k）;

创建辅助变量和函数。

ik = @（x）英国\（lk \ x）;ik_m = @（x）-ik（m * x）;瞬间= IK（F）;

定义一个与前三个矩相对应的频率序列。在这里,S.是频率。

Syms.S.锭队= [1 s ^ 2 s ^ 4];

设置第一时刻，与之相同d_Static在前面的计算中。

矩= d_static;

计算剩下的瞬间。

为了p=2:numel(powers) momentPre = iK_M(momentPre);为了pp = 1：磁幂（动量）矩（pp）=泰勒（动量次（pp），a，a，a，absply，'命令'，3）;结尾时刻= c * momentPre;时刻=(时刻;时刻);结尾

结合矩项来创建近似解析频率响应Happox．

HApprox =斯保尔*时刻;

显示前三个时刻。因为表达式很大，所以您可以仅部分显示它们。

momentString = arrayfun (@char时刻,'统一输出'、假);freqString = arrayfun (@char。’,'统一输出'、假);表(freqString momentString,'variablenames'，{'弗里克斯'那'时刻'}）

ans =3×2表FreqPowers时刻  _______________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________ {' 1 '}{”——(5 ^ (1/2)* (220 * 200 + * * cos (8352332796509007/9007199254740992) + 116 * 5 ^ (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函数一种和S.．您可以使用Gearting的MATLAB功能而无需符号数学工具箱。

Happoxfun = Matlabfunction（Happox，'vars', [s]);

比较纯粹的数字和象征性导出的频率响应

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

freq = 2 * pi * logspace（0,1）;

计算频率响应的数值a = afixed * perturbfactor，也就是说，对于周围的小扰动一种．

摄动因子= 1 + 0.25;KFixed =双(潜艇(K, AFixed * perturbFactor));MFixed =双(潜艇(M, AFixed * perturbFactor));H_Numeric = 0(大小(频率));为了1:numel(freq) sVal = 1i*freq(k);H_Numeric(k) = c*((sVal^2*MFixed + KFixed)\F);结尾

计算频率范围内频率响应的符号值a = perturbfactor * afix．

H_Symbolic = HApproxFun (AFixed * perturbFactor,频率* 1);

绘制结果。

图重对数(频率/(2 *π),abs (H_Symbolic),频率/(2 *π),abs (H_Numeric),'k *'）;Xlabel（'频率'）;ylabel（'频率响应'）;传奇(“象征”那'数字'）;

分析和数字曲线以所选择的间隔靠近，但通常，曲线在邻域之外的频率发散s = 0.．此外，对于价值观一种远离ablix，泰勒近似引入了更大的错误。