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

这个例子展示了如何找到参数化位移的解析表达式,联合的一个微不足道的悬臂桁架结构在静态和频率域。静态情况下精确的解析表达式。频率响应函数的表达式是一个近似降维版本的实际频率响应。

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

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 = 1 e2;Eval = 1 e7;

指定其他桁架钢筋的截面面积。

AFixed = 1;

接下来,定义本地特定桁架元素的刚度矩阵。

计算当地的刚度矩阵k并将其表示为一个象征性的函数。力和位移的桁架元素通过当地的刚度矩阵有关。桁架元素的每个节点有两个自由度,水平和垂直位移。每个桁架的位移矩阵对应元素是一个列[x_node1、y_node1 x_node2, y_node2]。由此产生的刚度矩阵是一个4×4矩阵。

信谊一个Elt真正的G = (cost罪(t) 0 0;0 0 cost罪(t)];kk = * 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 {因为 (t)}^{2}}{l} \\ σ_{3} = \frac{一个 E {罪 (t)}^{2}}{l} \end{array}$

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

信谊ρ6毫米= *ρ* l / * [2 1; 1 2];m =简化(G *毫米* G);localMassFn = symfun (m, [t、lρ]);

现在,定义桁架的酒吧作为一个矩阵酒吧。指数的酒吧、关节,开始和结束的关节是如下图所示。

矩阵的行数酒吧桁架是酒吧的数量。的列酒吧有四个元素,代表以下参数:

长度(l)
角度对水平轴(t)
指数的联合
指数结束联合

酒吧= 0 (2 * N + 2 * (N - 1), 4);

定义上和对角酒吧。注意,感兴趣的酒吧(N + 1)届酒吧或酒吧4号从左边第一斜杆。

为n = 1: n lelem = L / n;酒吧(n:) = (lelem 0 n, n + 1);lelem =√(L / N) ^ 2 + H ^ 2);酒吧(N + N:) = (lelem,量化(H、L / N)、N + 1 + N, N + 1);结束

定义低和垂直的酒吧。

为n = 1: 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

系统的质量矩阵米和系统刚度矩阵K是象征性的矩阵。将从每个酒吧局部元素矩阵组装成相应的全球矩阵。

K =符号(0 (numDofs numDofs));M =符号(0 (numDofs numDofs));为j = 1:尺寸(酒吧、1)%结构刚度和质量矩阵对酒吧j。lelem =酒吧(j, 1);24 =酒吧(j, 2);kelem = localStiffnessFn (24、lelem Eval);melem = localMassFn (24、lelem rhoval);n1 =酒吧(j, 3);n2 =酒吧(j, 4);%的酒吧不是参数化区域内,设置刚度%和质量矩阵的数值。注意,尽管这些值% kelem和melem没有象征性的参数,它们仍然是%符号对象(象征性的数字)。如果j ~ = N + 1 kelem =潜艇(kelem, AFixed);melem =潜艇(melem AFixed);结束%安排指数。第九= (n1 * 2 - 1, n1 * 2, n2 * 2 - 1, n2 * 2);%将局部元素嵌入到系统全局矩阵。K(第九,ix) = K(第九,ix) + kelem;M(第九,ix) = M(第九,ix) + melem;结束

消除附加自由度在墙上。

wallDofs = [1、2、2 * (N + 1) + 1, 2 * (N + 1) + 2];%自由度来消除:K (wallDofs) = [];K (:, wallDofs) = [];米(wallDofs:) = [];米(:,wallDofs) = [];

F是负载的负载向量中指定负的Y方向上右边的关节。

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

提取Y位移在上面右边的联合,创建一个选择向量。

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

从精确的符号创建Simscape方程静态情况下的解决方案

找到并绘制位移的解析解d的函数一个。在这里,K \ F代表了关节和位移c选择特定的位移。首先,显示符号代表数值使用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]);持有在;包含(的横截面,);ylabel (的位移,d ');标题(”)

图包含一个坐标轴对象。坐标轴对象包含截面,ylabel位移,d包含一个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。

AParamValues = AFixed * (1:5)”;d_NumericArray = 0(大小(AParamValues));为k = 1:元素个数(AParamValues) beginDim = (k - 1) *大小(k, 1) + 1;endDim = k *大小(k, 1);%创建一个数字刚度矩阵和提取数字解决方案。KNumeric =双(潜艇(K, AParamValues (K)));d_NumericArray (k) = c * (KNumeric \ F);结束

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

d_SymArray =双(潜艇(d_Static, AParamValues));

想象一个表中的数据。

T =表(AParamValues d_SymArray d_NumericArray)

T =5×3表1 AParamValues d_SymArray d_NumericArray其* * * -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)。
找到Pade逼近的传递函数 $H (年代, 一个) = c ({年代}^{2} 米 (一个) + K (一个))^{- - - - - - 1} F$ 在频率s = 0使用前三个系列的时刻。我们的想法是,给定一个传递函数,可以计算Pade逼近的时刻 $c (- - - - - - K (一个)^{- - - - - - 1} 米 (一个))^{j} K (一个)^{- - - - - - 1} F$ ,在那里 $j \in {0, 2, 4, 6, 。。。}$ 对应的幂级数 ${1, {年代}^{2}, {年代}^{4}, {年代}^{6}, 。。。}$ 。对于这个示例,计算HApprox(年代)前三项之和。这种方法是一个非常简单的技术,减少了传递函数的顺序。
进一步加快计算,近似每一刻的内部术语术语的泰勒级数展开的一个周围AFixed。

使用vpa加快计算速度。

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

计算的LU分解K。

(路、英国)= lu (K);

创建辅助变量和函数。

英国本土知识= @ (x) \(路\ x);iK_M = @ (x)动力学(M * x);momentPre =翼(F);

定义一系列频率对应于前三个时刻。在这里,年代是频率。

信谊年代sPowers = (1 s ^ 2 s ^ 4);

设置第一时刻,这是一样的d_Static在前面的计算。

时刻= d_Static;

计算剩余的时刻。

为p = 2:元素个数(sPowers) momentPre = iK_M (momentPre);为页= 1:元素个数(momentPre) momentPre (pp) =泰勒(momentPre (pp), AFixed,“秩序”3);结束时刻= c * momentPre;时刻=(时刻;时刻);结束

结合目前条件创建近似分析频率响应Happrox。

HApprox =斯保尔*时刻;

显示前三个时刻。由于表达式是大,你可以显示他们只是部分。

momentString = arrayfun (@char时刻,“UniformOutput”、假);freqString = arrayfun (@char。’,“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 * - 0.000000000000015007141035503798552656353179406 * (- 1) ^ 2 - 0.000000000000045855285883001825767717087472451 '}

频率响应转换为包含数值的MATLAB函数一个和年代。您可以使用生成的MATLAB函数没有符号数学工具箱。

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

比较纯粹的数字和象征性地派生的频率响应

不同的频率0来1在logspace,然后将范围转换为弧度。

频率= 2 *π* logspace (0,1);

计算出频率响应的数值= AFixed * perturbFactor,小扰动一个。

perturbFactor = 1 + 0.25;KFixed =双(潜艇(K, AFixed * perturbFactor));MFixed =双(潜艇(M, AFixed * perturbFactor));H_Numeric = 0(大小(频率));为k = 1:元素个数(频率)sVal = 1我*频率(k);H_Numeric (k) = c * ((sVal ^ 2 * MFixed + KFixed) \ F);结束

计算符号值的频率范围和频率响应= perturbFactor * AFixed。

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

策划的结果。

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

图包含一个坐标轴对象。坐标轴对象包含频率、ylabel频率响应包含2线类型的对象。一个或多个行显示的值只使用这些对象标记代表符号,数字。

分析和数值曲线选择间隔很近,但是,总的来说,小区外的频率曲线偏离s = 0。同时,的值一个远离AFixed,泰勒近似引入了更大的错误。