第六章-有限元三角形单元程序设计
有限元用三角形单元分析优秀课件
引入支 承条件
解方程 求位移
求应力
现结合图 5-21 说明。图中离散化后的三角形网格，有 4 个单元和 6 个结点， 进行了编号。支承在结点 4，5，6 处有四个支杆支撑。
(1) 建立整体刚度矩阵
结构中有六个结点，共有 12 个结点位移分量和 12 个结点力分量。由结点位移 求结点力F：
F K
分块形式
(4-28)
平面三角形单元整体分析
2) 求各单元的贡献矩阵 K e
以单元②为例，贡献矩阵K 2 由式(4-41)求出：
Str 1 2 3 4
1
2
10
1 1
02
0 2
3
K 2 Et
4
4
1 0 1 2
31 13
5
00
2 0
10
1 1
6
1 2
5 6
01
1
00
2 1
2
②子块“搬家”——把 k e 中按局部码排列的 9 个子块“搬家”，变为 K e 中按总码排列的 9
个子块。 K e 中的其余 27 个子块，则用零子块来填充。
以单元②为例，局部码 1，2，4 对应于总码 2，4，5，按对应关系搬家，得出单元②的贡献矩阵 K 2 ：
总码 1 2
K 2 3
4 5 6
0 2
3 1
1 2 1 3 0 1
0 0 2 0 2 0
1 0 1 1 0 1
（t——板厚）
(4-26)
4.4 平面三角形单元分析
弹性力学平面问题三角形单元的刚度矩阵 k e 具有对称性和奇异性等特点。 三角形单元刚度矩阵 k e 可写成分块形式：
FF12e
有限元单元法程序设计
有限元单元法程序设计是有限元分析（FEA）中的重要环节，它通过将连续的物理问题离散化为大量的、相互之间仅按特定方式相互联系的有限个单元的组合，从而进行求解。
以下是一个简单的有限元单元法程序设计的例子：
1.定义节点和单元：首先，我们需要定义模型的节点（nodes）和单元（elements）。
节点是空间中的点，而单元是由节点连接而成的物理实体。
2.建立网格：然后，我们需要根据模型的形状和大小，建立起一个合适的网格。
这个网格应该能够捕捉到模型的主要特征，并且足够细以捕捉到细节。
3.定义材料属性：接下来，我们需要为每个单元定义材料属性，比如弹性模量、泊松比、密度等。
4.施加载荷和约束：然后，我们需要根据问题的要求，对模型施加载荷和约束。
例如，我们可能需要施加压力、重力等载荷，以及位移、转动等约束。
5.进行有限元分析：最后，我们使用有限元方法进行求解。
这包括计算每个节点的位移和应力，以及根据这些结果进行后处理，比如生成报告、生成可视化图像等。
以上就是一个简单的有限元单元法程序设计的过程。
在实际应用
中，还需要考虑很多其他的因素，比如模型的复杂性、计算资源的限制等。
因此，编写一个有效的有限元程序需要深入理解有限元方法、计算机科学和工程知识。
matlab有限元三角形单元编程
matlab有限元三角形单元编程
在MATLAB中进行有限元分析，可以使用其自带的有限元分析工具箱（FEATool）进行编程。
以下是一个简单的例子，演示如何使用三角形单元进行有限元分析：
1. 打开MATLAB，进入FEATool环境。
2. 创建新的有限元模型。
选择“Model”菜单下的“New Model”选项，进入“Model Builder”界面。
3. 在“Model Builder”界面中，选择“2D Triangle”单元类型，并在绘图区域中绘制出三角形网格。
4. 在“Model Builder”界面中，设置材料属性、边界条件和载荷等参数。
5. 运行有限元分析。
选择“Model”菜单下的“Solve”选项，进行有限元求解。
6. 查看结果。
选择“Model”菜单下的“Results”选项，可以查看位移、应力、应变等结果。
以上是一个简单的例子，演示了如何使用三角形单元进行有限元分析。
在实际应用中，还需要根据具体问题进行详细的建模和计算。
第6章——常应变三角形单元
[S1
S 2 S3 ]
ci 1-µ bi 2
bi E µb = Si DB = i i 2 A(1 − µ 2 ) 1-µ ci 2
µ ci
平面应变：用平面应变弹性矩阵代入得到类似结果。
22:42
有限单元法
崔向阳
16
单元应变和应力矩阵
由于同一单元中的D、B矩阵都是常数矩阵，所以S矩阵也是常 数矩阵。也就是说，三角形三节点单元内的应力分量也是常 量。 当然，相邻单元的E， µ， A和bi、ci(i，j，m)一般不完全相同， 因而具有不同的应力，这就造成在相邻单元的公共边上存在 着应力突变现象。但是随着网格的细分，这种突变将会迅速 减小。
m j
h
1 i U = ∫∫ (σ xε x + σ yε y + τ xyγ xy )hdxdy x 2 A 1 T T T T = ∫∫ σ T εhdxdy σ = ( D ε ) = ε D 2 A
−1
u1 u2 u 3
u ( x, y ) = {1 x
1 x1 y} 1 x2 1 x3
y1 y2 y3
−1
u1 u2 u 3
4
22:42
有限单元法
崔向阳
平面三角形单元
假设
{ N1
N2
N 3 } = {1 x
m 相邻单元的位移在公共边上是连续的 形函数在单元上的面积分和边界上的线积分公式为 i
j p
A ∫ ∫A Ni dxdy = 3
式中 lij 为 ij 边的长度。
Ni =1 i m j
22:42
1 ∫ij Ni dl = 2 lij
有限元单元法程序设计
有限元单元法程序设计有限元单元法是一种用于工程结构分析和设计的计算方法，它将大型结构分解为许多小的离散单元，通过分析单元之间的相互作用来预测结构的力学行为。
有限元单元法程序设计是指针对特定工程问题，编写计算机程序来实现有限元分析的过程。
下面将介绍有限元单元法程序设计的基本流程和关键要点。
一、问题建模和网格划分有限元单元法程序设计的第一步是对工程结构进行合理的建模和网格划分。
建模的目的是将实际结构抽象为适用于有限元分析的数学模型，包括定义结构的几何特征、材料属性、边界条件等。
网格划分是将结构分解为许多小的单元，每个单元具有一定的形状和尺寸，以便于数值计算。
常用的单元形状包括三角形、四边形、四面体、六面体等，根据结构的特点选择合适的单元形状和尺寸。
二、单元刚度矩阵和载荷矩阵的求解在有限元单元法程序设计中，需要编写算法来求解每个单元的刚度矩阵和载荷矩阵。
单元刚度矩阵描述了单元内部的力学性能，包括刚度、弹性模量、泊松比等，它们通常通过数学公式或有限元理论推导得到。
载荷矩阵描述了单元受到的外部荷载，可以是均匀分布载荷、集中载荷或者边界条件引起的约束力。
通过合适的数值积分方法，可以计算得到每个单元的刚度矩阵和载荷矩阵。
三、组装全局刚度矩阵和载荷向量在有限元单元法程序设计中，需要将所有单元的刚度矩阵和载荷向量组装成整个结构的全局刚度矩阵和载荷向量。
这涉及到单元之间的连接关系以及边界条件的处理。
采用适当的组装算法，可以将各个单元的刚度矩阵和载荷向量叠加在一起，形成整个结构的刚度矩阵和载荷向量。
四、求解位移和应力有限元单元法程序设计的最后一步是求解结构的位移和应力。
通过斯蒂芬-泰勒算法或者其他迭代算法，可以得到整个结构的位移分布，然后根据位移场计算各个点的应变和应力。
这一过程涉及到对整个结构刚度矩阵的求解和对位移的后处理。
有限元单元法程序设计是一个复杂而又精密的工作，需要深入理解有限元原理、结构力学知识和数学方法。
matlab有限元三角形单元程序
matlab有限元三角形单元程序以下是一个简单的 MATLAB 有限元三角形单元程序的示例：```matlab% 定义模型参数E = 1000; % 弹性模量nu = 0.3; % 泊松比thickness = 1; % 板的厚度% 定义节点坐标nodes = [0, 0; 1, 0; 1, 1; 0, 1; 0.5, 0.5];% 定义单元连接关系connectivity = [1, 2, 5; 2, 3, 5; 3, 4, 5; 4, 1, 5];% 计算总节点数和总单元数numNodes = size(nodes, 1);numElements = size(connectivity, 1);% 初始化全局刚度矩阵和载荷向量K = zeros(numNodes);F = zeros(numNodes, 1);% 循环遍历每个单元for i = 1:numElements% 查找当前单元的节点编号nodesIndex = connectivity(i, :);% 根据节点编号从全局坐标矩阵中取出节点坐标coordinates = nodes(nodesIndex, :);% 计算当前单元的局部刚度矩阵localK = calculateLocalStiffness(E, nu, thickness, coordinates);% 组装局部刚度矩阵到全局刚度矩阵中K(nodesIndex, nodesIndex) = K(nodesIndex, nodesIndex) + localK;% 计算当前单元的局部载荷向量localF = calculateLocalLoad(thickness, coordinates);% 组装局部载荷向量到全局载荷向量中F(nodesIndex) = F(nodesIndex) + localF;end% 边界条件：节点1固定K(1, :) = 0;K(1, 1) = 1;F(1) = 0;% 解线性方程组U = K \ F;% 输出位移结果disp('节点位移：');disp(U);% 计算应力结果stress = calculateStress(E, nu, thickness, nodes, connectivity, U);% 输出应力结果disp('节点应力：');disp(stress);% 计算局部刚度矩阵的函数function localK = calculateLocalStiffness(E, nu, thickness, coordinates)% 计算单元的雅可比矩阵J = (1/2) * [coordinates(2,1)-coordinates(1,1), coordinates(3,1)-coordinates(1,1);coordinates(2,2)-coordinates(1,2), coordinates(3,2)-coordinates(1,2)];% 计算雅可比矩阵的逆矩阵invJ = inv(J);% 计算单元刚度矩阵B = invJ * [-1, 1, 0; -1, 0, 1];D = (E/(1-nu^2)) * [1, nu, 0; nu, 1, 0; 0, 0, (1-nu)/2]; localK = thickness * abs(det(J)) * (B' * D * B);end% 计算局部载荷向量的函数function localF = calculateLocalLoad(thickness, coordinates) localF = zeros(3, 1);midPoint = [sum(coordinates(:,1))/3,sum(coordinates(:,2))/3];localF(3) = thickness * 1 * det([coordinates(1,:); coordinates(2,:); midPoint]);end% 计算各节点应力的函数function stress = calculateStress(E, nu, thickness, nodes, connectivity, U)stress = zeros(size(nodes, 1), 3);for i = 1:size(connectivity, 1)nodesIndex = connectivity(i, :);coordinates = nodes(nodesIndex, :);Ke = calculateLocalStiffness(E, nu, thickness, coordinates);Ue = U(nodesIndex);stress(nodesIndex, :) = stress(nodesIndex, :) + (Ke * Ue)';endstress = stress / thickness;end```这个程序实现了一个简单的平面三角形单元有限元分析，包括定义节点坐标和单元连接关系、计算全局刚度矩阵和载荷向量、施加边界条件、解线性方程组、计算节点位移和应力等。
2012硕研(平面问题三角形单元有限元程序设计)
0、
主程序段—1、2
REAL KS(200,100) （定义数组总刚度矩阵）
OPEN(5,FILE=‘FEMT3.IN’)(打开原始数据文件) DIMENSION LND(50,3),X(100),Y(100),JR(20,3),PJ(20,3),P(200)
OPEN(6,FILE=‘FEMT3.OUT’,STATUS=‘NEW’)（输出文件）
READ(5,*)((PJ(I,J),J=1,3),I=1,NPJ)（荷载信息。依 次为：荷载所在结点号，x、y方向荷载大小）
100
NW=0（半带宽） DO 110 IE=1,NE DO 110 I=1,3 DO 110 J=1,3 IW=IABS(LND(IE,I)-LND(IE,J)) IF(NW.LT.IW) THEN NW=IW ENDIF 110 CONTINUE NW=(NW+1)*2 IF(IPS.NE.0) THEN E=E/(1.0-PR*PR) PR=PR/(1.0-PR) ENDIF END
• 有关带宽，储存方式请自修《有限元法概论》 中§5-12；《结构矩阵分析原理》第五章。
6．子程序6（形成总刚度矩阵[KS]）
SUBROUTINE STIFF(NJ,NE,NJ2,NW,LND,X,Y,E,PR,T,KS) DIMENSION LND(NE,3),X(NJ),Y(NJ) REAL KS(NJ2,NW),KE(6,6) DO 5 I=1,NJ2 DO 5 J=1,NW 5 KS(I,J)=0. DO 10 IE=1,NE CALL STE(IE,NJ,NE,LND,X,Y,E,PR,T,KE)
半带宽计算示例：
半带宽d的一般公式 半带宽d=(相邻结点码的最大差值+1)×结点位移未知量数 三角形单元：半带宽d=(相邻结点码的最大差值+1)×2 图中相邻结点码的最大差值是4,故d=(4+1)×2=10
第六章三角形单元的有限元法程序设计
从第j列的主对角线元素起到该列上方第一个非零 元素为止，所含元素的个数称为第j列的列高，记为hj ； 如果把第j列上方第1个非零元素的行号记为mj，则第j 列的列高为 hj = j - mj + 1 其实，hj就是第j行的左带宽，因而必有 UBW= max(hj)
j=1,2, …,N
利用节点位移信息数组 ID （去约束后节点位移自 由度编码），可容易地确定刚度矩阵 [K] 任何一列的 列高。
{ } [S ]{ }
例1：对角受压的正方形薄板，载 荷沿厚度均匀分布，为2N/m。由 于对称性，取1/4部分作为计算对
2N/m y x
象，试用有限元程序进行计算。
2N/m 2m 2m
h 1.0m, E 1.0, 0.0
例2：简支梁，梁高3m,跨度18m，厚度1m，承 受均布荷载10N/m2。已知 E 2 1010 N / m2 , 0.167
y为挤压应力，属次要应 力。材料力学不考虑 , FEM与弹性力学误差较大。
6.2 提高计算精度的方法
（1） 计算结果的整理 计算结果包括位移和应力两个方面。在位移方 面，一般无须进行整理工作。应力结果则需要 整理。通常认为计算出的应力是三角形单元形 心处的应力。而相邻单元之间的应力存在突变， 甚至正、负符号都不相同。为了由计算结果推 算出结构内某一点的接接实际的应力，必须通 过某种平均计算。通常可采用两单元平均法或 绕结点平均法。
平均法整理单元应力
两单元平均法：把两个相邻单元中的常应力加 以平均，用来表示公共边界中点处的应力。 绕结点平均法：把环绕某一结点的各单元常应 力加以平均，用以表示该结点的应力。在内结 点效果较好，而在边界结点可能很差，一般改 为应由内结点的应力外推计算出来。 （2）网格的细分 通过网格的细分，使每个单元的面积缩小，那 么尽管每个单元是应变、常应力单元，仍可较 好地反映结构中的应力变化，使得到的解答收 敛于问题的精确解。
弹性力学有限元-06-三角形单元
第6章 有限元法
6.1 概述 The finite element method
在17世纪， 牛顿和莱布尼茨发明了积分法，证明了该运算具有 整体对局部的可加性。
在18世纪,著名数学家高斯提出了加权余值法及线性代数方程组 的解法。另一位数学家Lagrange提出泛函分析。泛函分析是将 偏微分方程改写为积分表达式的另一途经。
1960年，美国加州大学伯克利分校的R.W.Clough教授在论文 中提出了“有限单元”，这样的名词。 值得骄傲的是我国南京 大学冯康教授在此前后独立地在论文中提出了“有限单元”
第6章 有限元法
The finite element method
有限元计算方法作为一种技术更多的与FEM软件的发展紧密的 结合起来。方法不断更新，优胜劣汰，传承和发展.有限元算法 （FFE）。在传统有限元分析的数值计算方法之中，有直接计算 法（DirectSolver）与迭代法（Iterative 所谓快速解法）两种。
3、单元内任意点的位移列阵f
{ f } [u ]T
（6-3）
y
m
v
·
u
j
i
y
m
· j
i
x
x
图6-1
4、单元内任意点的应变列阵
{} [ x y xy ]T
（6-4）
5、单元内任意点的应力列阵
{} [ x y xy ]T
y
m
· j
i
（6-5）
6、几何方程（Geometric equax tions）
3
1
3
5
7
9
66 6
④
⑤
555
单元、节点需编号
8
三角形常应变单元程序的编制与使用共18页文档
三角形常应变单元程序的编制与使用有限元法是求解微分方程边值问题的一种通用数值方法，该方法是一种基于变分法（或变分里兹法）而发展起来的求解微分方程的数值计算方法，以计算机为手段，采用分片近似，进而逼近整体的研究思想求解物理问题。
有限元分析的基本步骤可归纳为三大步：结构离散、单元分析和整体分析。
对于平面问题，结构离散常用的网格形状有三角形、矩形、任意四边形，以三个顶点为节点的三角形单元是最简单的平面单元，它较矩形或四边形对曲边边界有更好Array的适应性，而矩形或四边形单元较三节点三角形有更高的计算精度。
Matlab语言是进行矩阵运算的强大工具，因此，用Matlab语言编写有限元中平面问题的程序有优越性。
本章将详细介绍如何利用Matlab语言编制三角形常应变单元的计算程序，程序流程图见图1。
有限元法中三节点三角形分析结构的步骤如下：1）整理原始数据，如材料性质、荷载条件、约束条件等，离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。
2）单元分析，建立单元刚度矩阵。
3）整体分析，建立总刚矩阵。
4）建立整体结构的等效节点荷载和总荷载矩阵5）边界条件处理。
6）解方程，求出节点位移。
7）求出各单元的单元应力。
8）计算结果整理。
计算结果整理包括位移和应力两个方面；位移计算结果一般不需要特别的处理，利用计算出的节点位移分量，就可画出结构任意方向的位移云图；而应力解的误差表现在单元内部不满足平衡方程，图1 程序流程图单元与单元边界处应力一般不连续，在边界上应力解一般与力的边界条件不相符合。
1.1 程序说明% 三角形常应变单元求解结构主程序●功能：运用有限元法中三角形常应变单元解平面问题的计算主程序。
●基本思想：单元结点按右手法则顺序编号。
●荷载类型：可计算结点荷载。
●说明：主程序的作用是通过赋值语句、读取和写入文件、函数调用等完成算法的全过程，即实现程序流程图的程序表达。
1 程序准备format short e %设定输出类型clear all %清除所有已定义变量clc %清屏●说明：format short e －设定计算过程中显示在屏幕上的数字类型为短格式、科学计数法；clear all －清除所有已定义变量，目的是在本程序的运行过程中，不会发生变量名相同等可能使计算出错的情况；clc －清屏，使屏幕在本程序运行开始时2 全局变量定义global NNODE NPION NELEM NVFIX NFORCE COORD LNODS YOUNG POISS THICKglobal FORCE FIXEDglobal BMATX DMATX SMATX AREAglobal ASTIF ASLOD ASDISPglobal FP1说明：NNODE—单元结点数，NPION—总结点数， NELEM—单元数，NVFIX—受约束边界点数，NFORCE—结点力数，COORD—结构结点坐标数组，LNODS —单元定义数组，YOUNG—弹性模量，POISS—泊松比，THICK—厚度FORCE —节点力数组(n,3) n:受力节点数目,(n,1):作用点,(n,2):x方向,(n,3):y方向； FIXED—约束信息数组(n,3) n:受约束节点数目, (n,1):约束点 (n,2)与(n,3)分别为约束点x方向和y方向的约束情况,受约束为1否则为0BMATX—单元应变矩阵(3*6)， DMATX—单元弹性矩阵(3*3)，SMATX—单元应力矩阵(3*6)，AREA—单元面积ASTIF—总体刚度矩阵，ASLOD—总体荷载向量，ASDISP—结点位移向量FP1—数据文件指针3 打开文件FP1=fopen('input.txt','rt'); %打开输入数据文件存放初始数据●说明：FP1=fopen('input.txt','rt'); －打开已存在的输入数据文件input.txt，且设置其为只读格式，使程序在执行过程中不能改变输入文件中的数值，并用文件句柄FP1来执行FP2=fopen('output.txt','wt'); －打开输出数据文件，该文件不存在时，通过此命令创建新文件，该文件存在时则将原有内容全部删除。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
1
0
1
1
0
0
0
12
2(1)
*
材料矩阵计算程序
4%
((1)*(1-2*))* ... [(1) 0 0 0; (1) 0 0 0;
(1) 0 0 0; 0 0 0 (1-2*)/2 0 0; 0 0 0 0 (1-2*)/2 0; 0 0 0 0 0 (1-2*)/2];
G
E 2(1
)
,
(x123)/3; (y123)/3;
()
三角形单元程序
Ue 12deTKede
*
三角形单元程序
%
1 (i)((i));
0.5*(x1*y22*y33*y11*y32*y13*y2); %
2*2;
(12)*[(y23) (y31) (y12)]; % .
(12)*[(x32) (x13) (x21)]; % .
B(0, 2)
C(4, 1)
*
三角形单元程序
;
4; 2;
1.0;
0.0;
's
1;
%
*
A(0, 0) B(0, 2)
D(4,0) F 1 C(4, 1)
16; 8; 2**; 3; 2; (1)*(1); *; *;
% % %
三角形单元程序
x0=[]; 11
11 x0=[x0; (1)*
-0.5**(1+(1))*(1-(1))];
*
材料矩阵计算程序
D
E
1 2
1
0
1 0
0
0
1
2
材料矩阵计算程序
2% ((1)*(1-2*))* ... [(1) 0;
(1) 0; 0 0 (1-2*)/2];
3% ((1)*(1-2*))* ... [(1) 0; (1) 0; (1) 0; 0 0 0 (1-2*)/2];
1
D E(1) (1)(12)
B矩阵计算程序
*
[]2()
%: % %: %% k% -.
= (); 1
(i); 1
(j); ()()();
刚度矩阵组装程序
*
边界条件施加程序（置一法）
function [kk,ff]=feaplyc2(kk,ff,bcdof,bcval)
%---------------------------------------------------------% Purpose: % Apply constraints to matrix equation [kk]{x}={ff} % Variable Description: % kk - system matrix before applying constraints % ff - system vector before applying constraints % bcdof - a vector containging constrained d.o.f % bcval - a vector containing contained value %-----------------------------------------------------------
*
(,1); (); (); (,1); (,1); (,3); (,3); (,1); (3); (3,3);
(1);
*
三角形单元程序
三角形单元程序
1
(1)(,1); (2)(,2); (3)(,3);
%1 () %2 () %3 ()
x10((1),1); y10((1),2); 1 x20((2),1); y20((2),2); 2 x30((3),1); y30((3),2); 3
Hale Waihona Puke 1E(1
2
)
2G
0 0 0
2G
0
0
0
D
对
2G 0 0 0
G
0
0
称
G 0
G
*
局部单元系统自由度计算程序
[]()
%: % %: % - "“ %%-
0; 1 = ((i)-1)*; 1 1; (k);
*
[2]2d()
%: % % %: %%- x %- y
1 i1=(1)*2+1; i21+1; 2(11)(i); 2(22)(i); 2(31)(i); 2(32)(i);
三角形单元程序设计
问题描述
考虑一个平面应力问题如图所示，假设厚度1，材料为各项同 性，杨氏模量为1，泊松比为ν=0，相关力和位移边界条件如 图中所示，问题左端为固定约束。试用两个三角形单元分析 此问题，三角形单元的网格划分如图所示。试求问题各节点 位移u、v和应力σx，σy和σ。
A(0, 0)
F 1 D (4, 0)
();
*
三角形单元程序
Ni
1 2A
ai
bi x
ci
y
0.5*(x1*y22*y33*y11*y32*y13*y2); 2*2;
Ni,x
(12)*[(y23) (y31) (y12)]; x (12)*[(x32) (x13) (x21)]; y
Ni,y
22d(); 2'**2*;
Ke BTDBhdxdy BTDBA A
(3); %以下是求各个节点的应力值，如果节点被几个单元共有，则取这些单元中应力的平 均值
1 (); 1 {}(i); 1:3 ()()();
()();
*
三角形单元程序
% ('01','w'); (,'" "\n'); (,'"x" "y" "u" "v" "" "" ""\n'); (,' "", %88, \n'); 1
*
三角形单元程序
() > . > []; 1
1 [; (1)*(1) (1)* (1)*(1)1;]; [; (1)* (1)*1 (1)*(1)1;];
*
三角形单元程序
[]; []; 11
[ 1+2*(1) 2+2*(1)]; [ 0 0];
A(0, 0) B(0, 2)
D (4, 0) C(4, 1)
22d(); %
2*;
%
*;
%
1:3
()(i); %
()(i); %
0.5*'***;
*
ε(x,y)B (x,y)de
σ(x,y)D ε(x,y)
U 12ijijd
= (,1); (1); 1
1:3 (())(())+1; {()}((()));
三角形单元程序
%系统中每个编号共用次数 % 用来存储每个编号节点出现的单元号构成矩阵来构成所有的节点情况
2();
*
1; % (,1); % []2();
%
\; [ ](); \; \;
*
三角形单元程序
K ~df
dK~1f
0.5*'*1*;
%
0;
1
%
(1)(,1); % 1 ()
(2)(,2); % 2 ()
(3)(,3); % 3 ()
x10((1),1); y10((1),2) 1 x20((2),1); y20((2),2) 2 x30((3),1); y30((3),2) 3
(,'%16.616.616.616.616.616.616.6e\n'0(i,1)0(i,2), (2*1)(2*i)(1)(2)(3));
1 (,'%888d\n'(i,1)(i,2)(i,3));
*
[]()
%: % %: %% - 's % 1% 2% 3% 4-
1% (1*)* ... [1 0; ... 1 0; ... 0 0 (1)/2];