Simscape悬臂桁架结构分析模型

打开生活的脚本

该算例说明了如何从静力和频域求出普通悬臂桁架结构节点位移的参数化解析表达式。静态情况的解析表达式是精确的。频率响应函数的表达式是实际频率响应的近似降阶版本。

这个示例使用了以下Symbolic Math Toolbox™功能:

使用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 = [cos(t) sin(t) 0 0;0 0 cost 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 {因为 （ t ）}^{2}}{l} \\ σ_{3.} ＝ \frac{一个 E {罪 （ t ）}^{2}}{l} \end{array}$

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

信谊ρmm = A*rho*l/6*[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 leelem = 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 leelem = 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);对于不在参数化区域A内的杆，设置刚度%和质量矩阵的数值。注意，尽管值%勒姆和肉搏战没有符号参数，他们仍然是%符号对象(符号数字)。如果j ~= N+1 kelem = subs(kelem,A，固定);melem =潜艇(melem AFixed);结束%排列索引。第九= (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];要消除的dfs百分比: 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 ');标题(＇＇）

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

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

信谊dssEq = simscapeEquation (d, d_Static);

显示for表达式的前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:numel(AParamValues) beginindim = (k-1)*size(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表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。

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

创建辅助变量和函数。

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

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

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

设定第一个时刻，和d_Static在前面的计算中。

时刻= d_Static;

计算剩下的力矩。

为p=2:numel(sPowers) momentPre = iK_M(momentPre);为pp=1:numel(momentPre) = taylor(momentPre(pp)，A, fixed，“秩序”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 + * * 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 '}*(a - 1)^2 - 0.000000000000045855285883001825767717087472451'}

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

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

比较纯数字和符号推导的频率响应

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

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

计算频率响应的数值= AFixed * perturbFactor，即对周围的一个小扰动一个。

微扰因子= 1 + 0.25;KFixed =双(潜艇(K, AFixed * perturbFactor));MFixed =双(潜艇(M, AFixed * perturbFactor));H_Numeric = 0(大小(频率));为sVal = 1i*freq(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 (的频率响应);传奇(“象征”，“数字”);

图中包含一个坐标轴。轴线包含2个线型对象。这些对象代表符号，数字。

在选定的区间内，解析曲线和数值曲线是相近的，但一般来说，曲线在邻域以外的频率上是发散的s = 0。的值一个远离AFixed，泰勒近似引入了更大的误差。

Simscape悬臂桁架结构分析模型

定义模型参数

将棒组装成符号全局矩阵

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

比较数值静态解和符号静态解万博尤文图斯

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

比较纯数字和符号推导的频率响应

符号数学工具箱文档

万博1manbetx

数学建模与符号数学工具箱

Simscape悬臂桁架结构分析模型

定义模型参数

将棒组装成符号全局矩阵

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

比较数值静态解和符号静态解万博 尤文图斯

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

比较纯数字和符号推导的频率响应

符号数学工具箱文档

万博1manbetx

数学建模与符号数学工具箱

比较数值静态解和符号静态解万博尤文图斯