有限元 第二次作业

2-2 图示悬臂板,属于平面应力问题,其网格图及单元、节点编号见图2-1,E=2、1×1011,u=0、28,演算其单刚阵到总刚阵得组集过程,并用MATLAB 软件计算总刚阵。

图2-1

答:根据图2-1所示列出单元节点列表:

i j k 1 3 5 4 2 2 5 3 3 2 6 5 4

1

6

2

(1)计算单元刚度阵 单元1得刚度矩阵: ,; 单元2得刚度矩阵:,; 单元3得刚度矩阵:,; 单元4得刚度矩阵:,; 总刚度矩阵:

[]???????

??

?

?????

?????++++++++++++=4

6,636,635

,642

,632,641

,636

,535

,525,515,514,523

,513,53

2,522,515,414,413,425

,315,314,323

,313,322

,346,236,235

,225,223

,242

,232,222,241,246,142

,141

,10

00000

00

000k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k K 节点

单元

Matlab 程序语言得编写:

function Idex

global gNode gElement gMaterial

gNode=[0、0 0、01

0、5 0、01

1、0 0、01

1、0 0、0

0、5 0、0

0、0 0、0]

%gNode 同样就是一个矩阵,每一行表示一个结点,第1 列就是结点得x 坐标,第2 列就是结点得y坐标

gElement=[3 4 5

2 3 5

2 5 6

1 2 6 ];

%gElement 就是一个矩阵,每一行表示一个单元,第1 行就是单元得第1 个结点号,第2 行就是单元得第2个结点号。

Return

function k=StiffnessMatrix(ie)

%计算单元刚度矩阵函数

global gNode gElement

k=zeros(6,6); %6x6单元刚阵

E=2、1*10^11; %材料特性

u=0、28 ; %材料特性

t=0、01; %材料特性

xi=gNode(gElement(ie,1),1);

yi=gNode(gElement(ie,1),2);

xj=gNode(gElement(ie,2),1);

yj=gNode(gElement(ie,2),2);

xm=gNode(gElement(ie,3),1);

ym=gNode(gElement(ie,3),2); %计算节点坐标分量

ai=xj*ym-xm*yj;

aj=xm*yi-xi*ym;

am=xi*yj-xj*yi;

bi=yj-ym;

bj=ym-yi;

bm=yi-yj;

ci=-(xj-xm);

cj=-(xm-xi);

cm=-(xi-xj);

d=[1,xi,yi;1,xj,yj;1,xm,ym];

area=det(d); %计算单元面积

B=[bi 0 bj 0 bm 0 ;0 ci 0 cj 0 cm;ci bi cj bj cm bm];

B=B/2/area;

D=[1 u 0;u 1 0;0 0 (1-u)/2];

D=D*E/(1-u^2);

k=transpose(B)*D*B*t*abs(area); %计算单元刚度矩阵

Return

function gK=AssembleStiffnessMatrix

% 计算总刚阵

global gElement gK ie

gK=zeros(12,12);

for ie =1:1:4 %单元循环

k=StiffnessMatrix(ie);

for i=1:1:3 %节点循环

for j=1:1:3 %节点循环

for p=1:1:2 %自由度循环

for q=1:1:2 %自由度循环

m=(i-1)*2+p; %每个节点有2个自由度,i节点得第p个自由度为(i-1)*2+p

n=(j-1)*2+q; %每个节点有2个自由度,i节点得第p个自由度为(i-1)*2+p

M=(gElement(ie,i)-1)*2+p;

N=(gElement(ie,j)-1)*2+q;

gK(M,N)=gK(M,N)+k(m,n);

end

end

end

end

end

Return

则单元1得刚度矩阵为

>> StiffnessMatrix(1)

ans =

1、0e+010 *

2、0508 0 -2、0508 0、0410 0 -0、0410

0 5、6966 0、0319 -5、6966 -0、0319 0

-2、0508 0、0319 2、0531 -0、0729 -0、0023 0、0410

0、0410 -5、6966 -0、0729 5、6974 0、0319 -0、0008

0 -0、0319 -0、0023 0、0319 0、0023 0

-0、0410 0 0、0410 -0、0008 0 0、0008 单元2得刚度矩阵

>> StiffnessMatrix(2)

ans =

1、0e+010 *

2、0531 -0、0729 -2、0508 0、0319 -0、0023 0、0410

-0、0729 5、6974 0、0410 -5、6966 0、0319 -0、0008 -2、0508 0、0410 2、0508 0 0 -0、0410

0、0319 -5、6966 0 5、6966 -0、0319 0

-0、0023 0、0319 0 -0、0319 0、0023 0

0、0410 -0、0008 -0、0410 0 0 0、0008 单元3得刚度矩阵为

>> StiffnessMatrix(3)

ans =

1、0e+010 *

0、0023 0 -0、0023 0、0319 0 -0、0319

0 0、0008 0、0410 -0、0008 -0、0410 0

-0、0023 0、0410 2、0531 -0、0729 -2、0508 0、0319

0、0319 -0、0008 -0、0729 5、6974 0、0410 -5、6966

0 -0、0410 -2、0508 0、0410 2、0508 0

-0、0319 0 0、0319 -5、6966 0 5、6966

单元4得刚度矩阵

>> StiffnessMatrix(4)

ans =

1、0e+010 *

2、0531 -0、0729 -2、0508 0、0319 -0、0023 0、0410

-0、0729 5、6974 0、0410 -5、6966 0、0319 -0、0008

-2、0508 0、0410 2、0508 0 0 -0、0410

0、0319 -5、6966 0 5、6966 -0、0319 0

-0、0023 0、0319 0 -0、0319 0、0023 0

0、0410 -0、0008 -0、0410 0 0 0、0008

总刚度矩阵为

ans =

1、0e+011 *

Columns 1 through 8

0、2053 -0、0073 -0、0002 0、0041 0 0

0 0

-0、0073 0、5697 0、0032 -0、0001 0 0

0 0

-0、0002 0、0032 0、4106 -0、0073 -0、0002 0、0041 0 0

0、0041 -0、0001 -0、0073 1、1395 0、0032 -0、0001

0 0

0 0 -0、0002 0、0032 0、2053 0 -0、2051 0、0041

0 0 0、0041 -0、0001 0 0、5697 0、0032 -0、5697

0 0 0 0 -0、2051 0、0032 0、

2053 -0、0073

0 0 0 0 0、0041 -0、5697 -0、0073 0、5697

0 0 -0、4102 0、0073 0 -0、0073 -0、0002 0、0032

0 0 0、0073 -1、1393 -0、0073 0 0、0041 -0、0001

-0、2051 0、0041 0 -0、0073 0 0

0 0

0、0032 -0、5697 -0、0073 0 0 0

0 0

Columns 9 through 12

0 0 -0、2051 0、0032

0 0 0、0041 -0、5697

-0、4102 0、0073 0 -0、0073

0、0073 -1、1393 -0、0073 0

0 -0、0073 0 0

-0、0073 0 0 0

-0、0002 0、0041 0 0

0、0032 -0、0001 0 0

0、4106 -0、0073 -0、0002 0、0041

-0、0073 1、1395 0、0032 -0、0001

-0、0002 0、0032 0、2053 0

0、0041 -0、0001 0 0、5697

2-3 在平面问题有限元分析中,

(1)用到了哪些弹性力学中得基本方程？

答:平衡微分方程、几何方程、相容方程(形变协调方程)。

(2)力得平衡条件就是如何满足得？

答:根据能量守恒原理,有外力所作虚功应该等于内力虚功。也就就是结构在外载荷作用下处于平衡状态则在结构上得力在任意虚功位移上所作得虚功之与等于零。

以下就是用到得方程:

(3)变形协调条件就是如何满足得？

答:对材料进行线弹性与各向同性得假设,用弹性力学中应力-应变之间得关系得到变形协调条件。下面就是形变协调方程。

2-4 在平面三角形单元中得位移、应变、应力具有什么特征？

位移特征:(1)必须包含单元得刚体位移;

(2)必须包含单元得常应变状态;

(3)必须保证不偏惠各坐标轴;

(4)必须保证单元内位移连续。

应力特征:(1)三角形单元其应力仅与单元材料与几何尺寸有关,与节点位移有关,而与单元内位置坐标无关,也即这类单元内得应力就是常量。

(2)三角形单元内应力连续,但在公共边界上应力有突变,密布网格可以减少这种

冲突得不合理性。

应变特征:由于简单三角形单元取线性位移模式,其应变矩阵为常数矩阵,即在这样得位移模式下,三角形单元内得应变为某一常量。

2-5 在平面三角形单元中,当尺寸逐步缩小,单元中得位移、应变、应力具有什么特征？

当单元尺寸逐步减小时,单元各点得应变趋于相等,这时常量应变成为主要成分,因此,位移应能反应这种常应变状态,由于应力矩阵也就是常数矩阵,单元应力也就是常量。但就是相邻单元一般具有不同得力,在单元得公共边上会有应力突变,随着单元尺寸得逐步减小,这种突变会急剧降低,从而不会妨碍有线单元法得简答收敛于精确解。

清华大学有限元大作业

??????? ????? ???グ?? ? ??? ????φ??? ????φ??? ????φ??? ??? 1 Corresponding author: zzq_890709@https://www.wendangku.net/doc/2e914023.html,

?? Contents 0 ????????????? (2) 0.1 ?????? (2) 0.2 ????????? (2) 0.3 ?????? (3) 1????????仈 (4) 1.1 ???? (4) 1.2 ?????????????? (5) 1.3 ???? (7) 1.4 ???????ウ (13) 2 ???????????仈 (14) 2.1 ???? (14) 2.2 ???? (15) 2.3 ???? (16) 3 й???????????仈 (19) 3.1 ??????? (19) 3.2 ?????????????? (19) 3.3 ???? (20) 3.4 ???????ウ (22) 4 ??????????仈 (23) 4.1 ???? (23) 4.2 ???? (24) 4.3 ???? (25) 5 ??????????????仈 (32) 5.1 ???? (32) 5.2 ???? (33) 6 ???仈 (39) 6.1 ???? (39) 6.2 ??傼? (41) 6.3 ???? (44) 6.4 ???? (45) 7 ?? (48)

0 ????????????? 0.1 ?????? ????????????????лй??僔? 1????——?????????????Gambit. 2??????——????Matlab??????? 3????——??Matlab??.dat?????Tecplot??? 0.2 ????????? ?0 Matlab????????? ?????bandwidth.m ????????? elem_B_plane.m ????????B elem_D_plane.m ??????D elem_MK_plane.m ????????K elem_stressIntp.m ??儈??????? input.m ?????? input_readload.m ???? integrate_gauss_1.m 1?儈??? integrate_gauss_2.m 2?儈??? integrate_hammer.m й??Hammer?? main_AssemSolve.m ???????????Ka = P main_CalStress.m ?????????? output2tec.m ??.dat???Tecplot?? plot_mesh.m ?Matlab????? plot_stress.m ?Matlab??????? shape_cubic_820.m 3?8?20??Serendipity??????shape_lagrange.m 1?Lagrange???? shape_quadrangle_9.m ???9???????? shape_quadrangle_48.m ???4?8??Serendipity??????shape_triangle_36.m й??3?6????????

有限元上机实验指导书

弹性力学及有限元实验上机指导书 土木工程教学部 20１5年6月 一、ANSYS软件安装及有限元建模方法演示１、实验目的 掌握Ansys商用有限元软件的安装及了解有限元建模方法２、实验任务 １）Ansys商用软件8.1安装过程详解。

２）采用有限元直接建模法创建杆系模型演示。 3）采用有限元间接建模法创实体模型演示。 模型一：平面桁架如下图所示，长度单位为m，求支座反力和各杆内力。弹性模量2.1e+11，泊松比0.3，杆件截面面积均为0.01m2。 1/6 1/3 1/2 图1 平面桁架模型 模型二：正方形带孔平板，边长1m，小孔直径0.1m，板厚0.05m。弹性模量2.1e+11，泊松比0.3。上下边受均匀压力1000N。 图2 带孔正方形平板 3、实验方法 实验方法同课堂操作演示。 命令流(模型一) /PREP7 ET,1,LINK1 R,1,0.01, , MP,EX,1,,2.1e+11

MP,PRXY,1,,0.3 n,1,0,0 n,2,1,0 n,3,1,0.5 n,4,2,0 n,5,2,0.833 n,6,3,0 n,7,3,1 n,8,4,0 n,9,4,0.833 n,10,5,0 n,11,5,0.5 n,12,6,0 e,1,2 e,1,3 e,2,3 e,2,4 e,3,5 e,4,5 e,3,4 e,4,6 e,5,7 e,6,7 e,4,7 e,6,8 e,7,9 e,8,9 e,7,8 e,8,10

e,9,11 e,10,11 e,8,11 e,10,12 e,11,12 save 命令流(模型二) 二、利用ANSYS创建杆系或实体结构有限元模型 １、实验目的 掌握有限元建模的基本方法并能创建简单的杆系结构和实体结构有限元模型。 ２、实验任务 （1）采用直接建模法创建上节平面桁架结构模型。 （2）采用间接建模法创建上节带孔平板实体模型。 3、实验方法 实验方法同操作演示。 三、有限元求解及结果后处理演示 １、实验目的 掌握基本参数的设置、荷载施加方法及后处理操作。 ２、实验任务 （1）读入数据文件（命令流）的形式生成杆系结构有限元模型 （2）实常数、材料参数、求解参数设置演示 （3）位移约束、集中荷载施加方法演示 （4）计算求解与后处理操作演示。 3、实验方法 实验方法同课堂操作演示 附：后处理GUI操作及命令流操作 A、通过后处理提取节点计算结果 三种后处理操作： 1、plot results

重庆大学有限元第二次作业(刘静老师)

【有限元分析技术】第二次作业 科 目： 有限元分析技术 教 师： 姓 名： 学 号： 班 级： 类 别： 学术型 上课时间： 2016 年 11 月至 2017 年 1 月 考 生 成 绩： 卷面成绩 平时成绩 课程综合成绩 阅卷评语： 阅卷教师 (签名) 大学研究生院

第一章 题目概况 1.1 原始数据 矩形板尺寸如下图，板厚为5mm ，弹性模量为522.010/E N mm =? ，泊松比为0.27μ= 图1.1 原始计算简图 1.2工况选择 （1）试按下表的载荷约束组合，任选2种进行计算，并分析其位移、应力分布的异同。 表1 两种不同工况的载荷及约束 序号 载荷 约束 备注 1 向下均布载荷P=5N/mm,作用于ab 边 c ，d 点固定 2 向下均布载荷P=5N/mm,作用于ab 边 a ，b 点固定 3 向下均布载荷P=5N/mm,作用于ab 边 a ，c 边固定 还可讨论a ，c 点固定 4 向下均布载荷P=5N/mm,作用于cd 边 c ，d 点简支 5 向下均布载荷P=5N/mm,作用于cd 边 a ，b 点简支 6 向下均布载荷P=5N/mm,作用于cd 边 a ，c 边固定 还可讨论a ，c 点固定 7 向下集中载荷F=1000N,作用于ab 边中点 c ，d 点简支 8 向下集中载荷F=1000N,作用于ab 边中点 a ，b 点简支 9 向下集中载荷F=1000N,作用于ab 边中点 a ，c 边固定 还可讨论a ，c 点固定 10 向下集中载荷F=1000N,作用于cd 边中点 c ，d 点简支 11 向下集中载荷F=1000N,作用于cd 边中点 a ，b 点简支 12 向下集中载荷F=1000N,作用于cd 边中点 a ，c 边固定 还可讨论a ，c 点固定 1.3 工况选择结果及分析任务 （1）工况选择结果 根据表1的工况，选取工况1，2，8进行对比分析，选取结果如表2所示，为了方便下文中分别将序号1、2、8的工况称为工况一、工况二、工况三。 表2 分析工况的载荷及约束 序号 载荷 约束 备注 1 向下均布载荷P=5N/mm,作用于ab 边 c ，d 点固定 工况一 2 向下均布载荷P=5N/mm,作用于ab 边 a ，b 点固定 工况二 8 向下集中载荷F=1000N,作用于ab 边中点 a ， b 点简支 工况三

《结构分析中的有限元法》2015-有限元习题-参考答案

本科有限元习题参考答案

2015年3月10日作业 1、简述力学课程中介绍的各种力学模型的简化条件、基本假设和适用范围（包括有拉压杆模型、弯曲梁模型、平面应力和平面应变模型、轴对称模型、板模型、壳模型等） 2、给出弹性力学问题中平衡方程、几何方程、物理方程的表达式及其意义。 （1）平衡方程：

zy yz xz zx yx xy z yz xz z y xy zy y x zx yx x f y x z f x z y f z y x ττττττττσττσττσ====+??+??+??=+??+??+??=+??+??+??,000， 物理意义：应力分量与体力分量之间的关系。 （2）几何方程： z u x w y w z v x v y u z w y v x u zx yz xy z y x ??+??=??+??=??+??=??=??=??=γγγεεε,,,, 物理意义：应变分量与位移分量之间的关系。 （3）物理方程： [] [] [] zx zx yz yz xy xy y x z z z x y y z y x x G G G E E E τγτγτγσσμσεσσμσεσσμσε1,1,1) (1 ) (1 )(1 ===+-=+-=+-= 物理意义：应变分量与应力分量之间的关系。 3、简述最小势能原理的主要内容和主要公式。 根据虚功原理得到：??=-Γ T Ω T T 0Td Γδu d Ω）F δu -σδε（，由 )(21εδσεδδεU T T =?? ? ??=则0)21((=Γ-Ω-=∏??ΩΓ）Td u d F u T T T p σεδδ 其中，??ΩΓ Γ-Ω-=∏Td u d F u T T T p )21 (σε即为系统的总势能，它是弹性体变 形势能和外力势能之和。上面变分为零式表明：在所有区域内满足几何关系，在边界上满足给定位移条件的可能位移中，真实位移使系统的总势能取驻值(可证

华科大有限元分析题及大作业题答案——船海专业(DOC)

姓名：学号：班级：

有限元分析及应用作业报告 一、问题描述 图示无限长刚性地基上的三角形大坝，受齐顶的水压力作用，试用三节点常应变单元和六节点三角形单元对坝体进行有限元分析，并对以下几种计算方案进行比较： 1)分别采用相同单元数目的三节点常应变单元和六节点三角形单元计算； 2)分别采用不同数量的三节点常应变单元计算； 3)当选常应变三角单元时，分别采用不同划分方案计算。

二、几何建模与分析 图1-2力学模型 由于大坝长度>>横截面尺寸，且横截面沿长度方向保持不变，因此可将大坝看作无限长的实体模型，满足平面应变问题的几何条件；对截面进行受力分析，作用于大坝上的载荷平行于横截面且沿纵向方向均匀分布，两端面不受力，满足平面应变问题的载荷条件。因此该问题属于平面应变问题，大坝所受的载荷为面载荷，分布情况及方向如图1-2所示，建立几何模型，进行求解。 假设大坝的材料为钢，则其材料参数：弹性模量E=2.1e11，泊松比σ=0.3 三、第1问的有限元建模 本题将分别采用相同单元数目的三节点常应变单元和六节点三角形单元计算。 1）设置计算类型：两者因几何条件和载荷条件均满足平面应变问题，故均取Preferences为Structural 2）选择单元类型：三节点常应变单元选择的类型是PLANE42（Quad 4node42），该单元属于是四节点单元类型，在网格划分时可以对节点数目控制使其蜕化为三节点单元；六节点三角形单元选择的类型是PLANE183（Quad 8node183），该单元属于是八节点单元类型，在网格划分时可以对节点数目控制使其蜕化为六节点单元。因研究的问题为平面应变问题，故对Element behavior（K3）设置为plane strain。 3）定义材料参数 4）生成几何模 a. 生成特征点 b.生成坝体截面 5）网格化分：划分网格时，拾取所有线段设定input NDIV 为10，选择网格划分方式为Tri+Mapped，最后得到200个单元。 6)模型施加约束： 约束采用的是对底面BC全约束。 大坝所受载荷形式为Pressure，作用在AB面上，分析时施加在L AB上，方向水平向右，载荷大小沿L AB由小到大均匀分布（见图1-2）。以B为坐标原点，BA方向为纵轴y，则沿着y方向的受力大小可表示为： ρ（1） = gh P- =ρ g = - 10 {* } 98000 98000 (Y ) y

大作业报告参考2有限元学习心得

有限元学习心得 吴清鸽车辆工程 50110802411 短短八周的有限元课已经结束。关于有限元，我一直停留在一个很模糊的概念。我知道这是一个各个领域都必须涉及的点，只要有关于CAE分析的，几乎都要涉及有限元。总体来说，这是一门非常重要又有点难度的课程。 有限元方法(finite element method) 或有限元分析(finite element analysis)，是 求取复杂微分方程近似解的一种非常有效的工具，是现代数字化科技的一种重要 基础性原理。将它用于在科学研究中，可成为探究物质客观规律的先进手段。将 它应用于工程技术中，可成为工程设计和分析的可靠工具。本课程教学基本内容 有固体力学和结构力学简介；有限元法基础；桁架、梁、刚架、二维固体、板和 壳、三维固体的有限元法；建模技术；热传导问题的有限元分析；PATRAN软件 的使用. 通过有限元分析课程学习使我了解和掌握了一些有限元知识： 1．简要了解二维和三维固体以及桁架、梁和板结构的三组基本力学方程，即表示位移-应变关系的几何方程，表示应力-应变关系的本构方程和表示内力-外力关系的平衡方程。 2．了解利用能量法形成有限元离散系统方程的基本原理，即哈密尔顿原理。掌握有限元分 析的基本方法及步骤，包括域的离散、位移插值、构造形函数、单元有限元方程 的建立、坐标变换、整体有限元方程的组装、整体有限元方程的求解技术。 3．具体深入的了解并掌握桁架结构、梁结构、刚架结构、二维固体、板和壳结构、三维固体的有限元法分析技术，包括他们具体的形函数构造，应变矩阵，局部坐标系和整体坐标系中的单元矩阵。各种结构的实例研究。 4．了解并掌握建立高质量建模所涉及的各种关键技术。包括单元类型的选择，单元畸形的限制，不同阶数单元混用时网格的协调性问题，对称性的应用（平面对称、轴对称、旋转对称、重复对称），由多点约束方程形成刚域及应用（模拟偏移、不同自由度单元的连接、网格协调性的施加）等，以及多点约束方程的求解。以PATRAN有限元通用软件为例了解一般商业有限元软件的组成及结构。掌握PATRAN软件的基本使用。利用PATRAN软件上机实践完成两个上机练习：刚架结构有限元分析和三维固体有限元分析。 课程的具体学习内容： 内容： 1.三节点三角形单元：单元分析、总刚度矩阵组装、引入约束条件修正总刚度 矩阵、载荷移置、方程求解； 2.四边形单元分析、四节点四面体单元分析、八节点六面体单元分析；

有限元 第二次作业教学提纲

2-2 图示悬臂板，属于平面应力问题，其网格图及单元、节点编号见图2-1，E=2.1×1011，u=0.28，演算其单刚阵到总刚阵的组集过程，并用MATLAB 软件计算总刚阵。 图2-1 答：根据图2-1所示列出单元节点列表： i j k 1 3 5 4 2 2 5 3 3 2 6 5 4 1 6 2 （1）计算单元刚度阵 单元1的刚度矩阵：[] ????? ?????=15,514 ,513 ,515,414,41 3 ,415,314,313,3 1k k k k k k k k k k ，[] ?????????? ??????????=000000 000 00000000 000000 14 ,514 ,513 ,515,414,41 3,41 5,31 4,31 3 ,31 k k k k k k k k k k ； 单元2的刚度矩阵：[] ??? ? ????? ?=25 ,523 ,522 ,525,323,322,325,223,222,22 k k k k k k k k k k ，[] ????????? ?????? ?? ???=00 000000000000000000 0024,523,522,525,323 ,322,32 5 ,22 3 ,22 2,22k k k k k k k k k k ； 节点 单元

单元3的刚度矩阵：[] ??? ? ????? ?=36 ,635 ,632 ,636,535,532,536,235,232,23 k k k k k k k k k k ，[] ????????? ? ????? ?? ???=36,635 ,632,636,535,532 ,536,23 5 .23 2,2300000000 00 0000000000000000k k k k k k k k k k ； 单元4的刚度矩阵：[] ???? ? ???? ?=46,642 ,641,646,242,241 ,246,142,141 ,14 k k k k k k k k k k ，[] ?? ? ??? ? ?? ? ??????????=46,641 ,646,242,241.246,142 ,141,140000 000000000000000000000000k k k k k k k k k ； 总刚度矩阵：[][][][][][]4 3 2 1 4 1 k k k k k K e e +++=∑== []??????? ?? ? ????? ?????++++++++++++=4 6,636,635 ,642 ,632,641 ,636 ,535 ,525,515,514,523 ,513,53 2,522,515,414,413,425 ,315,314,323 ,313,322 ,346,236,235 ,225,223 ,242 ,232,222,241,246,142 ,141 ,10 00000 00 000k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k K Matlab 程序语言的编写： function Idex global gNode gElement gMaterial gNode=[0.0 0.01 0.5 0.01 1.0 0.01 1.0 0.0 0.5 0.0

有限元作业第二次作业

土木工程专业 有限元第二次作业 姓名： 班级： 学号： 指导教师： 二〇一五年6月12日

习 题：平面应力问题的八节点等参元，已给定8个节点 的坐标。试查资料并论述： 1、单元中位移函数u （ξ，η），v （ξ，η）和单元节点位 移{δe }的关系式； 2、[ B ]矩阵的计算步骤和计算式； 3、单元刚度矩阵[ k e ]的一般计算方法和计算步骤； 4、论述相邻单元间公共边界上位移的连续性； 5、如果给定母单元中点A ， （ξ，η），怎样求实际单元中与 A ， 相对应的点A （x ，y ）；反之，如果给定实际单元中的点A （x ，y ），怎样求其在母单元中对应点A ， （ξ，η）？ 6、如果已经求解得到单元8个节点的位移值{δe }怎样求单 元中某一点B （x ，y ）的应力？ 实际单元 1 2 6 7 Y 1 2 43 67 8η= 1η=﹣1 母单元 ξ= 1 ξ=﹣1

解： 1、此题分两步进行： 单元位移场的表达： 如图1所示,在任意四边形的每边中间设一附加节点，则单元边界就变成二次曲线的了。如果直接在整体坐标系(),x y 下，像八节点矩形元那样，构造双二次多项式的位移插值函数，则因曲边四边形单元边界是二次曲线，故边界上的位移是()x y 或的五次多项式， 它不能由曲边上三个节点的位移分量唯一地决定，从而不能保证相邻两个单元在公共边上位移的协调条件，所以在整体坐标系(),x y 下构造完全协调的位移插值 函数是很困难的，利用坐标变换，可将曲边四边形单元变换成基本单元，如图2所示的在自然坐标(),ξη下具有边长为2的八节点正方形单元，自然坐标系(),ξη是外节点坐标值为±1的局部坐标系。在自然坐标系的单元上构造 协调的位移插值函数，其形状函数是较普通的，取位移分量为,ξη的双二次多项式, 即： 2222 12345678222 2910111213141516u a a a a a a a a v a a a a a a a a ξηξξηηξηξηξηξξηηξηξη?=+++++++??=+++++++?? （1-1） 利用8 个节点的16 个位移分量可唯一确定16 个待定常数1216,,a a a …,，若代入8个节点的局部坐标值，得： 图1：在总坐标系中具有二 次曲边的四边形单元 图2：在自然坐标系中的 曲边四边形的基本单元

有限元分析报告大作业

有限元分析》大作业基本要求： 1.以小组为单位完成有限元分析计算，并将计算结果上交； 2.以小组为单位撰写计算分析报告； 3.按下列模板格式完成分析报告； 4.计算结果要求提交电子版，一个算例对应一个文件夹，报告要求提交电子版和纸质版。 有限元分析》大作业 小组成 员： 储成峰李凡张晓东朱臻极高彬月 Job name ：banshou 完成日 期： 2016-11-22 一、问题描述 （要求：应结合图对问题进行详细描述，同时应清楚阐述所研究问题的受力状况 和约束情况。图应清楚、明晰，且有必要的尺寸数据。）如图所示，为一内六角螺栓扳手，其轴线形状和尺寸如图，横截面为一外 接圆半径为0.01m的正六边形，拧紧力F为600N，计算扳手拧紧时的应力分布 图1 扳手的几何结构 数学模型

要求：针对问题描述给出相应的数学模型，应包含示意图，示意图中应有必要的尺寸数据；

图 2 数学模型 如图二所示，扳手结构简单，直接按其结构进行有限元分析。 三、有限元建模 3.1 单元选择 要求：给出单元类型， 并结合图对单元类型进行必要阐述， 包括节点、自由度、 实常数等。） 图 3 单元类型 如进行了简化等处理，此处还应给出文字说

扳手截面为六边形，采用4 节点182单元，182 单元可用来对固体结构进行

二维建模。182单元可以当作一个平面单元，或者一个轴对称单元。它由4 个结点组成，每个结点有2 个自由度，分别在x,y 方向。 扳手为规则三维实体，选择8 节点185单元，它由8 个节点组成，每个节点有3 个自由度，分别在x，y，z 方向。 3.2 实常数 （要求：给出实常数的具体数值，如无需定义实常数，需明确指出对于本问题选择的单元类型，无需定义实常数。） 因为该单元类型无实常数，所以无需定义实常数 3.3材料模型 （要求：指出选择的材料模型，包括必要的参数数据。） 对于三维结构静力学，应力主要满足广义虎克定律，因此对应ANSYS中的线性，弹性，各项同性，弹性模量EX：2e11 Pa, 泊松比PRXY=0.3 3.4几何建模由于扳手结构比较简单，所以可以直接在ANSYS软件上直接建模，在ANSYS建 立正六 边形，再创立直线，面沿线挤出体，得到扳手几何模型 图4 几何建模

有限元 第二次作业

2-2 图示悬臂板,属于平面应力问题,其网格图及单元、节点编号见图2-1,E=2、1×1011,u=0、28,演算其单刚阵到总刚阵得组集过程,并用MATLAB 软件计算总刚阵。 图2-1 答:根据图2-1所示列出单元节点列表: i j k 1 3 5 4 2 2 5 3 3 2 6 5 4 1 6 2 (1)计算单元刚度阵 单元1得刚度矩阵: ,; 单元2得刚度矩阵:,; 单元3得刚度矩阵:,; 单元4得刚度矩阵:,; 总刚度矩阵: []??????? ?? ? ????? ?????++++++++++++=4 6,636,635 ,642 ,632,641 ,636 ,535 ,525,515,514,523 ,513,53 2,522,515,414,413,425 ,315,314,323 ,313,322 ,346,236,235 ,225,223 ,242 ,232,222,241,246,142 ,141 ,10 00000 00 000k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k K 节点 单元

Matlab 程序语言得编写: function Idex global gNode gElement gMaterial gNode=[0、0 0、01 0、5 0、01 1、0 0、01 1、0 0、0 0、5 0、0 0、0 0、0] %gNode 同样就是一个矩阵,每一行表示一个结点,第1 列就是结点得x 坐标,第2 列就是结点得y坐标 gElement=[3 4 5 2 3 5 2 5 6 1 2 6 ]; %gElement 就是一个矩阵,每一行表示一个单元,第1 行就是单元得第1 个结点号,第2 行就是单元得第2个结点号。 Return function k=StiffnessMatrix(ie) %计算单元刚度矩阵函数 global gNode gElement k=zeros(6,6); %6x6单元刚阵 E=2、1*10^11; %材料特性 u=0、28 ; %材料特性 t=0、01; %材料特性 xi=gNode(gElement(ie,1),1);

有限元08上机作业

调试书本26到30页程序 开列数组维数： DIMENSION LOC(4,3),CX(6),CY(6),IFIX(6),F(12), 1GK(12,12),STRES(4,3),BAK(4,3,6) 结点集中力输入： DO 10 I=1,ND 10 F(I)=0.0 F(2)=-1.0 数据文件输入： 6,4,12,6,1.0E0,0.0,1.0,0.0,1 3,1,2 5,2,4 2,5,3 6,3,5 0.0,2.0 0.0,1.0 1.0,1.0 0.0,0.0 1.0,0.0 2.0,0.0 1,3,7,8,10,12 数据文件输出： NN NE ND NFIX E ANU T GM NTYPE 6 4 12 6 0.1000E+01 0.000 1.000 0.0000E+00 1 NODE X-LOAD Y-LOAD 1 0.000000E+00 -0.100000E+01 2 0.000000E+00 0.000000E+00 3 0.000000E+00 0.000000E+00 4 0.000000E+00 0.000000E+00 5 0.000000E+00 0.000000E+00 6 0.000000E+00 0.000000E+00 NODE X-DISP Y-DISP 1 -0.879121E-15 -0.325275E+01 2 0.879121E-16 -0.125275E+01 3 -0.879121E-01 -0.373626E+00 4 0.117216E-1 5 -0.835165E-15 5 0.175824E+00 -0.293040E-15 6 0.175824E+00 0.263736E-15 ELEMENT X-STR Y-STR XY-STR 1 -0.879121E-01 -0.200000E+01 0.439560E+00

有限元法课后习题答案

1、有限元是近似求解一般连续场问题的数值方法 2、有限元法将连续的求解域离散为若干个子域,得到有限个单元,单元和单元之间用节点连接 3、直梁在外力的作用下,横截面的内力有剪力和弯矩两个. 4、平面刚架结构在外力的作用下,横截面上的内力有轴力、剪力、弯矩. 5、进行直梁有限元分析,平面刚架单元上每个节点的节点位移为挠度和转角 6、平面刚架有限元分析，节点位移有轴向位移、横向位移、转角。 7、在弹性和小变形下，节点力和节点位移关系是线性关系。 8、弹性力学问题的方程个数有15个，未知量个数有15个。 9、弹性力学平面问题方程个数有8，未知数8个。 10、几何方程是研究应变和位移之间关系的方程 11、物理方程是描述应力和应变关系的方程 12、平衡方程反映了应力和体力之间关系的 13、把经过物体内任意一点各个截面上的应力状况叫做一点的应力状态 14、9形函数在单元上节点上的值,具有本点为_1_.它点为零的性质,并且在三角形单元的任一节点上,三个行函数之和为_1_ 15、形函数是_三角形_单元内部坐标的_线性_函数,他反映了单元的_位移_状态 16、在进行节点编号时,同一单元的相邻节点的号码差尽量小. 17、三角形单元的位移模式为_线性位移模式_- 18、矩形单元的位移模式为__双线性位移模式_

19、在选择多项式位移模式的阶次时,要求_所选的位移模式应该与局部坐标系的方位无关的性质为几何_各向同性 20、单元刚度矩阵描述了_节点力_和_节点位移之间的关系 21、矩形单元边界上位移是连续变化的 1.诉述有限元法的定义 答： 有限元法是近似求解一般连续场问题的数值方法 2.有限元法的基本思想是什么 答： 首先，将表示结构的连续离散为若干个子域，单元之间通过其边界上的节点连接成组合体。其次，用每个单元内所假设的近似函数分片地表示求解域内待求的未知厂变量。 3.有限元法的分类和基本步骤有哪些 答： 分类： 位移法、力法、混合法；步骤： 结构的离散化，单元分析，单元集成，引入约束条件，求解线性方程组，得出节点位移。 4.有限元法有哪些优缺点 答： 优点：

有限元作业

有限元作业

有限元分析大作业 学院： 班级： 姓名： 学号： 日期：

试题一（对应第二章） 如图所示，有一受轴向拉伸载荷2000P N =作用的变截面杆件，在0x =处，杆件截面积为2020A mm =，在180x L mm ==处，杆件截面积为201102 A mm =，杆件弹性模量为200GPa ，泊松比为0.3，试建立该杆件的有限元模型，并计算端部位移。（在划分网格时，沿长度方向取三个等长度杆单元） x P 0A 012A L 解：计算分析 000()(1)2(1)2(1)2x x x x A P x A A x A L P x A L P x E EA L σσσε===-=-==-

[ ]00 000022()[ln(2)]ln 2ln(2)12x x x x P dx PL PL u x dx L x L L x x EA EA EA L ε===--=--??- ??? ?? 0() 1.3860.1242mm PL u L EA == 数学建模：将其用二维模型进行降维处理，分为四个节点，三个等长度单元。 后处理

读出最大应力：1.750*10^2mpa 则计算得到右端部位移u（L）=0.12683 轴向位移随杆长变化图如下：

试题二（对应第三章） 一正方形平板,尺寸为40 mm×40 mm,厚度为2 mm,板中央有直径为d的圆孔如下图所示,板材弹性模量为200GPa,泊松比为0.3,在板的左端和右端分别施加20 MPa的拉力载荷.试建立该平板的有限元模型，并分别计算圆孔直 d=5,10,15,20和25mm时，平板开孔应力集中系数。

ansys上机作业

实验一坝体的有限元建模及应力应变分析 一、实验目的： 1、掌握ANSYS软件基本的几何形体构造方法、网格划分方法、边界条件施加方法及各 种载荷施加方法。 2、熟悉有限元建模、求解及结果分析步骤和方法。 3、能利用ANSYS软件对结构进行有限元分析。 二、实验设备： 微机，ANSYS软件。 三、实验内容： 计算分析模型如图所示，分析坝体的应力、应变。 四、实验步骤： 1 进入ANSYS 程序→ANSYS →change the working directory into yours →input Initial jobname: dam 2设置计算类型 ANSYS Main Menu: Preferences →select Structural → OK 3选择单元类型

ANSYS Main Menu: Preprocessor →Element Type→Add/Edit/Delete →Add →select Solid Quad 4node 42 →OK (back to Element Types window)→Options… →select K3: Plane Strain →OK→Close (the Element Type window) 4定义材料参数 ANSYS Main Menu: Preprocessor →Material Props →Material Models →Structural →Linear→Elastic→Isotropic→input EX:2.1e11, PRXY:0.3→ OK 5生成几何模型 生成特征点 ANSYS Main Menu: Preprocessor →Modeling →Create →Keypoints →In Active CS →依次输入四个点的坐标：input:1(0,0),2(1,0),3(1,5),4(0.45,5)→OK 生成坝体截面 ANSYS Main Menu: Preprocessor →Modeling →Create →Areas →Arbitrary →Through KPS→依次连接四个特征点，1(0,0),2(1,0),3(1,5),4(0.45,5) →OK 生成坝体截面如图一 图一 6 网格划分 ANSYS Main Menu: Preprocessor →Meshing →Mesh Tool→(Size Controls) lines: Set →依次拾取两条横边:OK→input NDIV: 15 →Apply→依次拾取两条纵边:OK →input NDIV: 40 →OK →(back to the mesh tool window)Mesh: Areas, Shape: Quad, Mapped →

有限元答案

1.1有限单元法中“离散”的含义是什么？有限单元法是如何将具有无限自由度的连续介质问题转变成有限自由度问题的？位移有限元法的标准化程式是怎样的？ （1）离散的含义即将结构离散化，即用假想的线或面将连续体分割成数目有限的单元，并在其上设定有限个节点；用这些单元组成的单元集合体代替原来的连续体，而场函数的节点值将成为问题的基本未知量。 （2）给每个单元选择合适的位移函数或称位移模式来近似地表示单元内位移分布规律，即通过插值以单元节点位移表示单元内任意点的位移。因节点位移个数是有限的，故无限自由度问题被转变成了有限自由度问题。 （3）有限元法的标准化程式：结构或区域离散，单元分析，整体分析，数值求解。 1.3单元刚度矩阵和整体刚度矩阵各有哪些性质？各自的物理意义是什么?两者有何区别? 单元刚度矩阵的性质：对称性、奇异性（单元刚度矩阵的行列式为零）。 整体刚度矩阵的性质：对称性、奇异性、稀疏性。 单元Kij物理意义Kij即单元节点位移向量中第j个自由度发生单位位移而其他位移分量为零时，在第j个自由度方向引起的节点力。 整体刚度矩阵K中每一列元素的物理意义是：要迫使结构的某节点位移自由度发生单位位移，而其他节点位移都保持为零的变形状态，在所有个节点上需要施加的节点荷载。 2.2什么叫应变能？什么叫外力势能？试叙述势能变分原理和最小势能原理，并回答下述问题：势能变分原理代表什么控制方程和边界条件？其中附加了哪些条件？ (1)在外力作用下，物体内部将产生应力ζ和应变ε，外力所做的功将以变形能的形式储存起来，这种能量称为应变能。 (2)外力势能就是外力功的负值。 (3)势能变分原理可叙述如下：在所有满足边界条件的协调位移中，那些满足静力平衡条件的位移使物体势能泛函取驻值，即势能的变分为零 δΠp=δUε+δV=0 此即变分方程。对于线性弹性体，势能取最小值，即 δ2ΠP=δ2Uε+δ2V≧0 此时的势能变分原理就是著名的最小势能原理。 势能变分原理代表平衡方程、本构方程和应力边界条件，其中附加了几何方程和位移边界条件。 2.3什么是强形式？什么是弱形式？两者有何区别？建立弱形式的关键步骤是什么？ 等效积分形式通过分部积分，称式 ∫ΩC T(v)D(u)dΩ+∫ΓE T(v)F(u)dΓ为微分方程的弱形式，相对而言，定解问题的微分方程称为强形式。 区别：弱形式得不到解析解。 建立弱形式的关键步骤：对场函数要求较低阶的连续性。 2.4为了使计算结果能够收敛于精确解，位移函数需要满足哪些条件？为什么？ 只要位移函数满足两个基本要求，即完备性和协调性，计算结果便收敛于精确解。 2.6为什么采用变分法求解通常只能得到近似解？变分法的应用常遇到什么困难？Ritz法收敛的条件是什么？ （1）在Ritz 法中，N决定了试探函数的基本形态，待定参数使得场函数具有一定的任意性。如果真实场函数包含在试探函数之内，则变分法得到的解答是精确的；如果试探函数取自完全的函数序列，则当项数不断增加时，近似解将趋近于精确解。然而，通常情况下试探函数不会将真实场函数完全包含在内，实际计算时也不可能取无穷多项。因此，试探函数只能是真实场函数的近似。可见，变分法就是在某个假定的范围内找出最佳解答，近似性就源于此。 （2）采用变分法近似求解，要求在整个求解区域内预先给出满足边界条件的场函数。通常情况下这是不可能的，因而变分法的应用受到了限制。 （3）Ritz 法的收敛条件是要求试探函数具有完备性和连续性，也就是说，如果试探函数满足完备性和连续性的要求，当试探函数的项数趋近于无穷时，则Ritz 法的近似解将趋近于数学微分方程的精确解。 3.1构造单元形函数有哪些基本原则？形函数是定义于单元内坐标的连续函数。单元位移函数通常采用多项式，其中的待定常数应该与单元节点自由度数相等。为满足完备性要求，位移函数中必须包括常函数和一次式，即完全一次多项式。多项式的选取应由低阶到高阶，尽量选择完全多项式以提高单元的精度。若由于项数限制而不能选取完全多项式时，也应使完全多项式具有坐标的对称性，并且一个坐标方向的次数不应超过完全多项式的次数。有时为了使位移函数保持一定阶次的完全多项式，可在单元内部配置节点。然而，这种节点的存在将增加有限元格式和计算上的复杂性，除非不得已才加以采用。形函数应保证用它定义的位移函数满足收敛要求，即满足完备性要求和协调性条件。 3.1构造单元形函数有哪些基本原则？试采用构造单元的几何方法，构造T10 单元的形函数，并对其收敛性进行讨论。 通常单元位移函数采用多项式，其中的待定常数由节点位移参数确定，因此其个数应与单元节点自由度数相等。根据实体结构的几何方程，单元的应变是位移的一次导数。为了反映单元刚体位移和常应变即满足完备性要求，位移函数中必须包含常数项和一次项，即完全一次多项式。 3.3何谓面积坐标？其特点是什么？为什么称其为自然坐标或局部坐标？ （1）三角形单元中，任一点P(x,y)与其3个角点相连形成3个子三角形,其位置可以用下述称为面积坐标的三个比值来确定： L1=A1/A L2=A2/A L3=A3/A 其中A1,A2,A3分别为P23,P31,P12的面积。 （2）面积坐标的特点： a T3单元的形函数Ni就是面积坐标Li b面积坐标与三角形在整体坐标系中的位置无关。 c三个节点的面积坐标分别为节点1（1, 0, 0）、节点2（0, 1, 0）、节点3（0, 0, 1），形心的面积坐标为（1/3, 1/3, 1/3）。 d单元边界方程为Li=0(i=1,2,3) e在平行于23边的一条直线上，所有点都有相同的面积坐标L1（L1对应的三角形具有相同的高和底边），而且L1就等于此直线至23边的距离与节点1至23边的距离之比值。

有限元上机题

注：题中E表示材料弹性模量，μ表示泊松比，ρ表示密度。 一静力结构分析 1 如图1所示为普通订书钉，E=2.1×105MPa，μ=0.3，横截面尺寸为宽B=0.64mm，高H=0.402mm。当订书钉被压入纸张时，约需要120N的载荷，载荷均匀地分布在订书钉上部。以下面两种情况进行有限元分析。（单位：mm）（1）订入时A、B点为铰支条件； （2）订入时A、B点为固定约束。 图1 载荷和尺寸情况 2、小型铁路桥由横截面积为3250mm2的钢制杆件组装而成。一辆火车停在桥上，其载荷施加在桥梁两侧的桁架上，单侧的桁架如图2所示，等效载荷为F1，F2，E=2.1×105MPa，μ=0.3，ρ=7.8×103kg/m3。试计算位置R处由于载荷作用而沿水平方向移动的距离以及支反力，同时，分析各个节点的位移和非单元应力。 图2 铁路桥单侧桁架及载荷情况 3如图3所示，模型参数为：E=3.0×1010Pa，A1=30m，A2=10m，B=80m，t=20m，p=2200Pa。 有关风载的确定，按照海洋井架行业标准，有以下方法： 风压（Pa）=0.6115×风速（m/s）×高度系数×形状系数对于一般的海洋井架及建筑物，高度在30m左右，高度系数取为1.1，形状系数取为1.25，风速取为47.8m/s。换算出来后得到的风压为2200Pa。

图3 高层建筑物及其风载荷 4对于含裂纹体的结构及材料，若按照线弹性力学分析，会在裂纹的尖端处产生应力的奇异性，这时需要计算裂纹尖端处的应力强度因子（对于Ⅰ型裂纹，有K1=σ（πa）1/2），并以应力强度因子作为准则来对材料的裂纹是否扩展进行判断。图4所示为一块矩形平板，其边缘存在长为a的裂纹，板的两端承受拉应作用。利用结构上下的对成性，取矩形的一半建立有限元模型，完成看一力σ 下工作： （1）球裂纹的张角θ（在施加载荷前=0，θ=0） （2）沿直线AO，画出y方向应力σy沿x变化的曲线图。假设 σy= K1/(2πx)1/2利用回归方法估计K1。将计算结果与计算无限大平板的修正结果进行比较，其中基于无限大平板的修正关系式为K1=1.2σ0（πa）1/2）（3）在裂缝尖端处，进行网格的细化，重新求解（2）中的问题。 模型中的相关参数为 E=2.1×105MPa，μ=0.3，L=400mm，a=9.5mm，b=95mm，σ0=450MPa。 图4 承受拉升载荷的带裂纹平面结构

有限元第二次作业

2-2图示悬臂板，属于平面应力问题，其网格图及单元、节点编号见图2-1，E=×1011，u=，演算其单刚阵到总刚阵的组集过程，并用MATLAB软件计算总刚阵。 图2-1 答：根据图2-1所示列出单元节点列表： i j k 1354 2253 3265 4162 （1）计算单元刚度阵 单元1的刚度矩阵：[] ? ? ? ? ? ? ? ? ? ? = 1 5,5 1 4,5 1 3,5 1 5,4 1 4,4 1 3,4 1 5,3 1 4,3 1 3,3 1 k k k k k k k k k k，[] ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = 1 4,5 1 4,5 1 3,5 1 5,4 1 4,4 1 3,4 1 5,3 1 4,3 1 3,3 1 k k k k k k k k k k； 单元2的刚度矩阵：[] ? ? ? ? ? ? ? ? ? ? = 2 5,5 2 3,5 2 2,5 2 5,3 2 3,3 2 2,3 2 5,2 2 3,2 2 2,2 2 k k k k k k k k k k，[] ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = 2 4,5 2 3,5 2 2,5 2 5,3 2 3,3 2 2,3 2 5,2 2 3,2 2 2,2 2 k k k k k k k k k k； 节点 单元

单元3的刚度矩阵：[] ??? ? ????? ?=36 ,635 ,632 ,636,535,532,536,235,232,23 k k k k k k k k k k ，[] ????????? ? ????? ?? ???=36,635 ,632,636,535,532 ,536,23 5 .23 2,2300000000 00 0000000000000000k k k k k k k k k k ； 单元4的刚度矩阵：[] ???? ? ???? ?=46,642 ,641,646,242,241 ,246,142,141 ,14 k k k k k k k k k k ，[] ?? ? ??? ? ?? ? ??????????=46,641 ,646,242,241.246,142 ,141,140000 000000000000000000000000k k k k k k k k k ； 总刚度矩阵：[][][][][][]4 3 2 1 4 1 k k k k k K e e +++=∑== []??????? ?? ? ????? ?????++++++++++++=4 6,636,635 ,642 ,632,641 ,636 ,535 ,525,515,514,523 ,513,53 2,522,515,414,413,425 ,315,314,323 ,313,322 ,346,236,235 ,225,223 ,242 ,232,222,241,246,142 ,141 ,10 00000 00 000k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k k K Matlab 程序语言的编写： function Idex global gNode gElement gMaterial gNode=[