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

打开生活的脚本

此示例说明如何在静态和频域中找到普通悬臂桁架结构接头位移的参数化解析表达式。静态情况下的解析表达式是精确的。频率响应函数的表达式是实际频率的近似降阶版本应急响应。

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

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

定义模型参数

这个例子的目的是解析地联系位移D为均匀截面面积参数A.悬臂桁架结构中特定的杆。控制方程用 $M \frac{D^{2.} x}{D T^{2.}} + K x = F$ ．在这里,D由悬臂梁右上关节处的荷载产生。桁架在最左边的接合处附在墙上。

桁架杆由密度均匀的线弹性材料制成。横截面积A.以蓝色突出显示的条的长度（见图）可能会有所不同。所有其他参数，包括其他条的均匀横截面面积，都是固定的。通过使用小的线性位移假设来获得尖端的位移。

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

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

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

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

rhoval=1e2；Eval=1e7；

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

AFixed = 1;

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

计算当地的刚度矩阵K并将其表示为符号函数。桁架单元的力和位移通过局部刚度矩阵进行关联。桁架单元的每个节点有两个自由度，水平和垂直位移。每个桁架单元的位移是一个柱矩阵，对应于[x_节点1，y_节点1，x_节点2，y_节点2].得到的刚度矩阵为4乘4矩阵。

信谊A.ELT真正的G=[cos（t）sin（t）0；0cos（t）sin（t）]；kk=A*E/l*[1-1；-11]；k=simplify（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{A. E 罪 (2. T)}{2. L} \\ σ_{2.} = \frac{A. E {因为 (T)}^{2.}}{L} \\ σ_{3.} = \frac{A. E {罪 (T)}^{2.}}{L} \end{array}$

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

信谊rhomm = 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

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

K =符号(0 (numDofs numDofs));M =符号(0 (numDofs numDofs));对于j = 1:尺寸(酒吧、1)构造杆j的刚度和质量矩阵。lelem=钢筋（j，1）；telem=钢筋（j，2）；kelem=局部刚度fn（telem，lelem，Eval）；melem=局部质量fn（telem，lelem，rhoval）；n1=钢筋（j，3）；n2=钢筋（j，4）；对于不在参数化区域A内的杆，设置刚度%和质量矩阵的数值。注意，尽管值%勒姆和肉搏战没有符号参数，他们仍然是%符号对象(符号数字)。如果j~=N+1kelem=subs（kelem，A，AFixed）；melem=subs（melem，A，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（墙自由度，：）=[]；K（：，墙自由度）=[]；M（墙自由度，：）=[]；M（：，墙自由度）=[]；

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

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

提取Y置换在最右上关节处，创建一个选择向量。

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

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

找到并绘制位移的解析解D作为A.．在这里,K\F表示所有接头和接头处的位移C选择特定的位移。首先，展示用16位数字表示数值的符号解。

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

ans=
                 
                   $- \frac{0.00000001118033988749895 (504.7737197247586 {A.}^{2.} + 384.7213595499958 A. + 22.3606797749979)}{A. (1．28 A. + 0.8944271909999158)}$

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

图中包含一个轴对象。axis对象包含一个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(……”

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

的范围内比较精确解析解和数值解A.参数值。这些值形成一个序列AFixed来5*a固定以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 (,)可以有效地快速检查参数的影响A.而不必为每个新值进行可能昂贵的数值计算A.．然而，对于带有附加参数的大系统，通常不可能找到精确的封闭解。因此，您通常会近似此类系统的解。这个例子近似于零频率附近的参数频率响应如下:

使用可变精度算法加快计算速度(vpa).
求传递函数的Padé近似 $H (s, A.) = C (s^{2.} M (A.) + K (A.))^{- 1.} F$ 频率附近s = 0用这个系列的前三个时刻。其思想是，给定一个传递函数，你可以计算Padé近似矩为 $C (- K (A.)^{- 1.} M (A.))^{J} K (A.)^{- 1.} F$ ,在那里 $J \in {0, 2., 4., 6., ．．．}$ 对应于幂级数项 ${1., s^{2.}, s^{4.}, s^{6.}, ．．．}$ .对于本例，请计算HApprox(年代)这种方法是降低传递函数阶数的一种非常基本的方法。
为了进一步加快计算速度，请使用泰勒级数展开式近似每个力矩项的内项A.周围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);

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

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

设置第一个力矩，该力矩与d_Static在前面的计算中。

力矩=d_静态；

计算剩余力矩。

对于p=2:numel(sPowers) momentPre = iK_M(momentPre);对于pp=1:numel（momentPre）momentPre（pp）=泰勒（momentPre（pp），A，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 + * * 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函数A.和s。您可以在不使用符号数学工具箱的情况下使用生成的MATLAB函数。

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

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

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

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

计算频率响应的数值= AFixed * perturbFactor，也就是说，对于周围的小扰动A.．

微扰因子= 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 *’); xlabel(“频率”）;ylabel (的频率响应）;传奇(“象征”,“数字”）;

图中包含一个轴对象。轴对象包含两个类型为line的对象。这些对象代表符号，数字。

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