2节点3节点Timoshenko梁单元刚度阵质量阵matlab符号推导

pasuka

利用符号运算推导的2节点和3节点 Timoshenko梁单元刚度和质量矩阵，应该不存在剪切锁死的问题
公式推导什么的就不具体列举了，实在太基础了～

Matlab语言: 高亮代码由发芽网提供
01 % 2-node 3-node Timoshenko Beam formula
02 % version: 0.1
03 % author: T.Q.(Shanghai Huitao Info-Technology Co. Ltd.)
04 % date; 22/08/2012
05 % e-mail: pasuka at foxmail.com
06 % Ref: M. Petyt, Introduction to finite element vibration analysis, p115
07
08 % Copyright (c) 2012, Shanghai Huitao Info-Technology Co. Ltd.
09 % All rights reserved.
10 %
11 % Redistribution and use in source and binary forms, with or without
12 % modification, are permitted provided that the following conditions are
13 % met:
14 %
15 %     * Redistributions of source code must retain the above copyright
16 %       notice, this list of conditions and the following disclaimer.
17 %     * Redistributions in binary form must reproduce the above copyright
18 %       notice, this list of conditions and the following disclaimer in
19 %       the documentation and/or other materials provided with the distribution
20 %      
21 % THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
22 % AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23 % IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
24 % ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
25 % LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
26 % CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
27 % SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
28 % INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
29 % CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
30 % ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31 % POSSIBILITY OF SUCH DAMAGE.
32 clear all;clc;
33 syms D C L ro I A x real
34 % D  ---- bending rigidity
35 % C  ---- shear rigidity
36 % L  ---- length of beam
37 % ro ---- density
38 % I  ---- inertia
39 % A  ---- cross area
40 F = [1,x,x^2,x^3,x^4,x^5];
41 % 3-node Timoshenko Beam
42 % deflection
43 w = F-D/C*diff(F,x,2);
44 % rotation
45 phi = diff(F,x,1);
46 T3(1,:) = subs(w,x,0);
47 T3(2,:) = subs(phi,x,0);
48 T3(3,:) = subs(w,x,L/2);
49 T3(4,:) = subs(phi,x,L/2);
50 T3(5,:) = subs(w,x,L);
51 T3(6,:) = subs(phi,x,L);
52 % shape function
53 N3 = [w;phi/T3;
54 % bending
55 Bk3 = diff(phi,x,1)/T3;
56 % shear
57 Bs3 = (diff(w,x,1)-phi)/T3;
58 stif3 = D*int(transpose(Bk3)*Bk3,x,0,L)+C*int(transpose(Bs3)*Bs3,x,0,L);
59 mass3 = ro*int(transpose(N3)*diag([A,I])*N3,x,0,L);
60 % simple stiff matrix
61 % simplify(stif3)
62
63 % 2-node Timoshenko Beam
64 T2(1,:) = subs(w(1:4),x,0);
65 T2(2,:) = subs(phi(1:4),x,0);
66 T2(3,:) = subs(w(1:4),x,L);
67 T2(4,:) = subs(phi(1:4),x,L);
68 N2 = [w(1:4);phi(1:4/T2;
69 Bk2 = diff(phi(1:4),x,1)/T2;
70 Bs2 = (diff(w(1:4),x,1)-phi(1:4))/T2;
71 stif2 = D*int(transpose(Bk2)*Bk2,x,0,L)+C*int(transpose(Bs2)*Bs2,x,0,L);
72 mass2 = ro*int(transpose(N2)*diag([A,I])*N2,x,0,L);

tonnyw

为什么没有剪锁死的问题？好像你的积分还是全积分。

pasuka

tonnyw 发表于 2012-8-22 22:22
为什么没有剪锁死的问题？好像你的积分还是全积分。

梁的挠度项里面已经考虑了剪切的影响，附件是一种比较直观且简单的推导方式，结果和胡海昌院士以及Petyt书中的结果一致

tonnyw

Thanks, I need to digest it a little.