有限元单元刚度矩阵
有限元法与ANSYS技术-刚度矩阵
k
e
BiT BjT
D
Bi
Bj
Bm
t
kii k ji
kij k jj
kim k jm
(3-35)
BmT
kmi kmj kmm
其中 krs Br T DBs t
Et
4 1 2
bcrbrbss1122cbrrccss
br cs
1
2
cr bs
cr cs
1
2
brb s
Kij
Kim
Kin
K
K j1
K ji
K jj
K jm
K jn
考虑到[k]扩充以后，除了对应的i, j, m 双行和双列 上的九个子矩阵之外，其余元素均为零，故（3-33）式
中的单元位移列阵{}e2n×1 便可用整体的位移列阵 {}2n×1 来替代。这样，（3-33）式可改写为
k
Re
2n2n
2n1
2n1
把上式对N个单元进行求和叠加，得
N e1
e
ui
vi
u j
v j
um
T
vm
且假设单元内各点的虚位移为{f *}，并具有与真实位移 相同的位移模式。
故有
f N e
(c)
参照（3-13）式，单元内的虚应变{ *}为
B e
(d)
于是，作用在单元体上的外力在虚位移上所做的功可写
为
({ }e )T Re
(f)
而单元内的应力在虚应变上所做的功为
码的排序一致。各单元的节点力列阵经过这样的扩充之
后就可以进行相加，把全部单元的节点力列阵叠加在一
起，便可得到 (l)式所表示的弹性体的载荷列阵，即
单元刚度矩阵
单元刚度矩阵单元刚度矩阵是在结构力学中的一个重要概念，它是一个矩阵，用来表示刚性和结构特性。
它可以用来描述广泛的结构，如桥梁，大型建筑物和其他复杂的构造。
它的研究有助于更好地理解结构的运动和反应。
它也可以用来预测和控制结构的变形和损坏，从而减少结构建设过程中可能发生的各种问题。
单元刚度矩阵是一个n x n等阶矩阵，其中n是一个复杂结构中的单元数量。
它代表了单元之间的约束关系，表明它们如何互相影响。
这也就是所谓的单元刚度矩阵。
每个矩阵元素代表了任意两个单元之间的受拉或受压力的数量，可以用来计算结构中每一个单元之间的刚度和约束。
单元刚度矩阵有几种不同的类型，其中一种是静态刚度矩阵，它用来表示复杂结构在静态荷载作用下的刚度。
它可以用来预测荷载作用下结构变形的情况，并作出相应的改善。
另外一种是有限元分析，它可以用来对复杂结构在动态荷载下的变形，受力，反应，以及可能发生的结构破坏作出分析。
单元刚度矩阵的计算方法有很多。
有些是利用有限元分析的方法来进行的，也有些是直接从节点和单元的计算和配置来得出的。
有些方法只需简单地求解结构中一组特定问题，而另一些方法则要求对结构中所有部件进行复杂的数值计算。
单元刚度矩阵的计算可以帮助从两个角度来改善设计：一方面，单元刚度矩阵可以帮助改善结构运动的性能，另一方面，它可以帮助减少结构上可能发生的变形以及提高结构的耐久性。
单元刚度矩阵的计算和研究非常重要，现代的结构力学和建筑设计工程正在用这个技术来设计新型的可靠性更高，耐久性更强的建筑结构。
基于单元刚度矩阵的计算和研究，科学家们可以更好地理解结构力学，并减少建筑物的再建设和变形，以及可能发生的损坏。
总之，单元刚度矩阵的研究和计算存在着很多的优势。
现代的结构力学和建筑设计都需要用到它，以便更好地分析和控制结构的变形和损坏。
它的研究也有助于开发更安全，更高效的建筑结构，有助于结构力学中的其他方面的研究。
optistruct 单元刚度矩阵
optistruct 单元刚度矩阵
在有限元分析中，单元刚度矩阵（Element Stiffness Matrix）是用于描述一个单元对应的局部坐标系下的刚度性质的矩阵。
OptiStruct是一种常用的有限元分析软件，它也根据单元的几何形状和材料特性计算出单元的刚度矩阵。
单元刚度矩阵描述了单元受力和变形之间的关系，它可以用于计算整个结构的全局刚度矩阵。
OptiStruct使用几何非线性、材料非线性和接触等特性来计算单元刚度矩阵。
根据不同的单元类型（如线性、非线性、壳单元等），OptiStruct采用不同的方法和公式来计算单元刚度矩阵。
一般来说，单元刚度矩阵的计算需要考虑以下几个方面：
1. 几何刚度：单元的形状和尺寸对刚度矩阵的计算有影响，如线性单元的刚度矩阵与单元长度有关。
2. 材料性质：材料的弹性模量和泊松比等材料特性对刚度矩阵的计算有影响。
3. 边界条件：单元所在的整体结构的边界条件对刚度矩阵的计算也有影响。
4. 单元类型：不同的单元类型具有不同的刚度矩阵计算方法。
了解单元刚度矩阵的计算对于进行有限元分析模拟和结果预测非常重要。
通过OptiStruct等有限元分析软件，可以方便地计算出各种类型的单元刚度矩阵，并进一步分析结构的强度和刚度等性能。
有限元 第3讲补充_平面问题-整体刚度矩阵
整体刚度矩阵
通过以上组装过程可以得到组装整体刚度矩阵的一般规则： 1 ）结构中的等效节点力是相关单元结点力的叠加，整体 刚度矩阵的子矩阵是相关单元的单元刚度矩阵子矩阵的集成； 2）当整体刚度矩阵中的子矩阵K rs 中r=s时，该节点（节点r 或s）被哪几个单元所共有，则K rs 就是这几个单元的刚度矩阵 e 中的子矩阵 K rs 的相加。如 K 33 应该是单元①－④中对应子矩阵 (1) (2) (3) (4) 的集成，即 K33 K33 K33 K33 K33
0
0
0
0 1 0 2 (2) 0 3 K 0 4 0 5
式中： Fi (2) ——②号单元中第i（i=1,3,4）节点所受力；
K (2) ——②号单元的扩大刚度矩阵。
y
4 ④ ② ① 1 2 3③ 5
(1)
0 0
0 0
0 0 1 0 0 2 (1) 0 0 3 K 0 0 4 0 0 5
4 ④ ② ① 1 2 3③
5
x
o
K (1) ——①号单元的扩大刚度矩阵或称为单元贡献矩阵。
5
整体刚度矩阵
y
4 ④ ② ① 14 ④ ② ① 1 2 3③ 5
x
o
(1) (2) (1) (1) (2) (2) K11 K11 K12 K13 K13 K14 0 (1) (1) (3) (1) (3) (3) K 22 K 22 K 23 K 23 0 K 25 K 21 (1) (2) (1) (3) (1) (2) (3) (4) (2) (4) (3) (4) K 31 K 31 K 32 K 32 K 33 K 33 K 33 K 33 K 34 K 34 K 35 K 35 (2) (2) (4) (2) (4) (4) 0 K 43 K 43 K 44 K 44 K 45 K 41 (3) (3) (4) (4) (3) (4) 0 K K K K K K 52 53 53 54 55 55
汽车结构有限元分析03单元类型及单元分析
1.一维单元分析
主要有：杆单元、梁单元、管单元等 。
1.1杆单元---最简单的两节点一维单元， 用于杆件承受轴向力分析。
设杆单元横截面积为Ａ，长度为l，轴 向分布载荷q为(x) 。单元2个节点的位移 向量为： e ui u j T
由空间弹性力学几何方程，得应变表达式：
{} [B]{ }e [[B1 ][B2 ][B20 ]]{ }e
由空间弹性力学物理方程，单元内的应力可 以表示成：
[ ] [D][ ] [D][B]{ }e [S]{ }e
单元刚度矩阵为 ：
[k]e
[B]T [D][B]dV
[k1e1
[k
e 21
这其中设定单元位移模式，利用虚功 原理建立单元节点力与节点位移关系并组建 单元刚度矩阵的过程，我们将其称为单元分 析。
为了使有限元法的解在单元尺寸逐步趋 小时能够收敛于精确解，所构造的单元位移 函数必须满足以下三方面的条件：
1）位移模式中必须包括反映刚体位移的项；
2）位移模式中必须包括反映常应变的线性位 移项；
这样空间梁单元就由3个节点组成i,，j,k 点必
须在一个平面内，但不能共线。i节点到j节
点为单元坐标系的x轴，y轴(或z轴)在节点i、
j和k构成的平面上且与x轴垂直，应用右手定
则可以确定另一坐标iz, 轴j, k(或y轴)。
三点
确定后，单元坐标系即确定，梁单元的截面
方位也就完全确定下来。所增加的一个用于
] ]
[k1e2 ]
[k
e 22
]
[k1e20
[k
有限元法基础重点归纳(精)
30、有限元法的任务:建立和求解整个弹性体的节点位移和节点力之间的关系的平衡方程。31、单元刚度矩阵:表达了单元节点位移与节点力之间的转换关系。
32、单元刚度矩阵的性质:①单元刚度矩阵中每个元素有明确的物理意义②K e是对称矩阵③K e的每一行或每一列元素之和为零,因此K e为奇异矩阵④K e不随单元的平行移动或作n π角度的转动而改变。33、刚度集成法集成规律:①先对每个单元求出其单元刚度矩阵K e ,而且以分块形式按节点编号顺序排列②将单元刚度矩阵扩大阶数为2n*2n ,并将单元刚度矩阵中的子块按局部码与总码的对应关系,搬到扩大后的矩阵中,形成单元贡献矩阵K e。③将所有单元贡献矩阵同一位置上的分块矩阵简单叠加成总体刚度矩阵中的一个子矩阵,各行各列都按以上步骤即形成总体刚度矩阵K。34、整体刚度矩阵的性质:①整体刚度矩阵是对称矩阵②整体刚度矩阵中每一元素的物理意义:整体刚度矩阵的第一列元素代表使第一个节点在x方向有一单元位移,而其余节点位移皆为零时必须在节点上施加的里。对于K的其余各列也有类似意义③整体刚度矩阵K的主对角线上的元素总是正的④整体刚度矩阵K是一个稀疏阵⑤整体刚度矩阵K是一个奇异阵。35、带形矩阵:整体刚度矩阵K的非零元素分布在以主对角线为中心的斜带形区域内的矩阵。
γxy
=E 1−μ
2∗
1−μ2
γxy
42、制造位移函数:{u (x,y =α1+α2x +α3y
v (x,y =α4+α5x +α6y
43、等参单元精度比四边形单元高,四边形精度比三角形精度高。
44、轴对称问题:很多工程物件,它们的几何形状承受的载荷以及约束条件都对称于其一固定轴,这即为对称轴,此时载荷作用下的位移、应变和应力也对称于该对称轴的问题。45、等参数单元:优点:①形状方位任意,适应性好,精度高,容易构造高阶单元②具有统一形式,规律性强,采用数值积分算,程序处理方便③高阶等参单元精度高,描述复杂边界,形状能力强,所需单元少。缺点:①单元各方向尺寸要尽量接近②单元边界不能过于曲折,不能有拐点折点,尽量接近直线或抛物线③边之间夹角要尽量接近直角④单元形状不能过度畸变,边中节点不能过于偏离中间。46、有限元法基础理论:弹性力学,材料力学
c3d8有限元单元方程推导过程
有限元单元方程推导过程1.引言有限元分析是一种数值计算方法，用于求解结构力学、流体动力学等领域的物理问题。
在有限元分析中，有限元单元是构成整个有限元模型的基本单元，通过推导有限元单元的方程，可以实现对结构或系统的精确分析和计算。
本文将从有限元方法的基本原理出发，详细介绍有限元单元方程的推导过程。
2.有限元方法基本原理有限元方法是将连续的物理问题离散化，转化为有限个代表性元素的集合，通过对每个元素施加适当的边界条件和力学方程，最终得到整个系统的解。
有限元方法通过有限元单元之间的相互作用，从而模拟整个系统的行为。
3.有限元单元的概念有限元单元是有限元模型中最小的离散单元，它是对实际的结构或系统进行离散化的结果。
不同的物理问题和结构，可以采用不同类型的有限元单元进行离散化，如梁单元、壳单元、板单元等。
4.有限元单元方程的一般形式有限元单元方程的一般形式可以表示为：\[K_{e}U_{e}=F_{e}\]其中\(K_{e}\)为有限元单元的刚度矩阵，\(U_{e}\)为有限元单元的位移矢量，\(F_{e}\)为有限元单元的荷载矢量。
5.有限元单元方程推导的基本步骤有限元单元方程的推导主要包括以下几个基本步骤：5.1 单元刚度矩阵的推导首先需要根据有限元单元的几何形状和材料性质，推导出单元刚度矩阵。
单元刚度矩阵可以通过对单元内部的应变能量或者应力-应变关系进行积分得到。
5.2 单元位移矢量的表示在推导单元方程过程中，需要选择合适的位移矢量表示方式，可以采用基函数展开的方法，将位移矢量表示为一组未知系数乘以基函数的线性组合形式。
5.3 单元荷载矢量的求解单元荷载矢量是由外部施加的荷载和边界条件共同决定的，在推导单元方程的过程中需要将这些荷载转化为局部坐标系下的形式，并利用位移矢量的表示方式，将荷载矢量表达为位移矢量和未知系数的线性组合。
5.4 单元方程的组装需要将单元刚度矩阵、位移矢量和荷载矢量组装成完整的单元方程，可以通过坐标变换或者有限元单元之间的关系对单元方程进行组装。
有限元中单元刚度矩阵和残差矩阵关系(一)
有限元中单元刚度矩阵和残差矩阵关系(一)有限元中单元刚度矩阵和残差矩阵关系什么是有限元分析有限元分析是一种解决连续介质力学问题的数值计算方法。
它将复杂的结构划分为许多小的单元，并对每个单元进行离散化和建模，以求解分析问题。
在有限元分析中，单元刚度矩阵和残差矩阵是两个重要的概念。
单元刚度矩阵单元刚度矩阵是对每个单元在局部坐标系下进行局部建模后的刚度矩阵。
它描述了单元内部各个节点之间的刚度关系。
单元刚度矩阵的计算通常基于材料的性质、几何形状和边界条件等。
残差矩阵残差矩阵是在有限元分析中引入的一个重要概念，用于描述节点的约束关系。
它是根据边界条件和节点位移的计算结果生成的。
在有限元分析中，为了保证整个结构的连续性和平衡，必须对节点之间的位移进行限制和约束。
残差矩阵表示这些约束关系。
单元刚度矩阵和残差矩阵的关系单元刚度矩阵和残差矩阵之间存在着紧密的关系。
这种关系可以通过有限元分析的公式推导得到。
通常来说，在线性静力学问题的有限元求解过程中，可以通过以下步骤来求解整体的刚度矩阵和残差矩阵：1.将整个结构划分为若干个单元，对每个单元进行局部建模和刚度矩阵的计算。
2.根据节点之间的约束关系，将所有单元的局部刚度矩阵进行组装，得到整个结构的总体刚度矩阵。
3.在施加边界条件的情况下，求解整体刚度矩阵和施加边界条件的节点位移，得到节点位移的解。
4.根据节点位移的解，计算整体结构的残差矩阵，即节点受到的力的不平衡情况。
因此，可以说单元刚度矩阵是构建整体刚度矩阵的基础，而残差矩阵则是在整体刚度矩阵和节点位移的解的基础上得到的。
单元刚度矩阵和残差矩阵之间的关系是有限元分析中求解问题的关键所在。
以上就是有限元中单元刚度矩阵和残差矩阵关系的简述和解释。
通过对单元刚度矩阵和残差矩阵的理解和计算，可以帮助我们更好地理解和解决连续介质力学问题。
有限元法中空间梁单元矩阵的组集
有限元法中空间梁单元矩阵的组集有限元法是一种常用的工程分析方法，用于解决结构力学问题。
它将连续体分割成大量的小区域，即有限元，通过对每个小区域的应力和变形进行计算，进而得到整个结构体系的力学性能。
在有限元法中，梁单元是一种常用的元素类型，用于分析和设计梁结构。
空间梁单元是三维空间中的梁元素，用于分析和设计具有一定长度、截面和材料特性的结构。
它通常由两个节点和一些节点上的自由度组成。
在有限元法中，通过对空间梁单元的刚度矩阵进行组集，可以对梁结构进行强度和刚度等力学性能的计算和评估。
空间梁单元的刚度矩阵是一个重要的参数，它描述了梁元素在受力作用下的应力和变形关系。
在有限元法中，通过将梁单元划分为许多小单元，每个小单元的刚度矩阵可以通过材料性质和几何形状等参数进行计算和组集。
通过将所有小单元的刚度矩阵按照一定规则组合，可以得到整个梁单元的刚度矩阵。
组集空间梁单元的刚度矩阵需要考虑以下几个方面的因素：1. 几何因素：梁单元的几何形状对其刚度矩阵有很大的影响。
梁单元的长度和横截面形状将决定其自由度的数量和排列方式。
在组集刚度矩阵时，需要将这些几何因素考虑进去，确保计算结果的准确性。
2. 材料性质：梁单元的材料性质对其刚度矩阵的计算也有影响。
不同材料的弹性模量和剪切模量将导致不同的应力和变形响应。
在组集刚度矩阵时，需要将材料性质的影响考虑进去，以获得准确的结果。
3. 节点约束：梁单元的刚度矩阵还需要考虑节点约束的影响。
节点约束可以限制节点上的位移和旋转自由度，影响刚度矩阵的计算。
在组集刚度矩阵时，需要将节点约束的影响考虑进去，以获得准确的结果。
通过以上几个因素的考虑，可以得到空间梁单元的刚度矩阵。
该刚度矩阵可以用于计算梁单元的应力、变形、弯曲刚度、切割刚度等力学性能，也可以用于解决梁结构在受力作用下的静力分析和动力分析问题。
对于空间梁单元刚度矩阵的组集，可以按照如下步骤进行：1. 划分梁单元为若干个小单元：根据需要和几何形状，将梁单元划分为若干个小单元，每个小单元可以视为一个简单的线性结构，易于计算其刚度矩阵。
有限元单元刚度矩阵计算方法
有限元单元刚度矩阵计算方法
有限元单元刚度矩阵是有限元分析中的一个关键组成部分，它描述了结构中每个元素在承受载荷时的刚度响应。
以下是一个计算有限元单元刚度矩阵的基本步骤：
1. 确定元素类型和参数：首先需要确定所使用的元素类型（例如，杆、梁、板、壳等），以及这些元素的参数，如横截面面积、惯性矩、厚度等。
2. 建立局部坐标系：为每个元素建立一个局部坐标系。
在局部坐标系中，可以方便地描述元素内部的应力和应变。
3. 计算应变矩阵：根据有限元理论，计算元素两端的节点坐标差值，并由此得到应变矩阵。
4. 计算应力矩阵：根据材料的物理性质和胡克定律（Hooke's law），将应变矩阵转换为应力矩阵。
5. 形成刚度矩阵：将应力矩阵乘以相应的刚度系数，得到该元素的刚度矩阵。
6. 组装整体刚度矩阵：将所有元素的局部刚度矩阵组合起来，形成整体结构的刚度矩阵。
7. 施加边界条件和载荷：根据实际问题的边界条件和载荷，对整体刚度矩阵进行修正。
8. 求解线性方程组：通过求解修正后的线性方程组，得到结构中每个节点的位移。
以上步骤仅为有限元分析中的一种基本方法，实际应用中可能还需要考虑更多的因素，如非线性行为、材料失效等。
此外，有限元分析软件（如ANSYS、SolidWorks等）通常已经内置了这些计算过程，用户可以直接调用相应的功能进行有限元分析，而无需手动编写代码。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
实现功能：输入欲求的单元号,得到该单元的面积和该单元刚度矩阵；输入完所有的单元号,得到所有的单元刚度矩阵后,直接得到半带宽存储的数组。
（所有的单元刚度矩阵和半带宽数组分别存在所有单元的刚度矩阵.txt和SK矩阵.txt）第1单元面积为：0.5000弹性模量、泊松比和厚度分别为:100.0000 0.3000 0.1000单元1单元的应力矩阵-109.8901 -32.9670 109.8901 0.0000 0.0000 32.9670-32.9670 -109.8901 32.9670 0.0000 0.0000 109.8901-38.4615 -38.4615 0.0000 38.4615 38.4615 0.0000第1单元的单元刚度阵为7.4176 3.5714 -5.4945 -1.9231 -1.9231 -1.64843.5714 7.4176 -1.6484 -1.9231 -1.9231 -5.4945-5.4945 -1.6484 5.4945 0.0000 0.0000 1.6484-1.9231 -1.9231 0.0000 1.9231 1.9231 0.0000-1.9231 -1.9231 0.0000 1.9231 1.9231 0.0000-1.6484 -5.4945 1.6484 0.0000 0.0000 5.4945半带宽存储的数组7.4176 3.5714 -1.9231 -1.6484 -5.4945 -1.92317.4176 -1.9231 -5.4945 -1.6484 -1.9231 0.00009.8901 -2.3810 -4.9451 4.7619 -3.0220 -0.45799.8901 4.7619 -4.9451 -0.7326 0.5495 0.000020.8791 -2.3810 0.5495 -2.3810 -10.989 1.648413.7363 -2.3810 -3.0220 1.9231 -3.8462 0.00003.4341 1.1905 -0.9615 1.9231 0.0000 0.00005.2198 1.6484 -2.7473 0.0000 0.0000 0.000011.9506 -3.5714 0.0000 0.0000 0.0000 0.00006.5934 0.0000 0.0000 0.0000 0.0000 0.0000#include <stdio.h>#define jds 5#define d 2main(){int nj,ne,i,j,k,a,n,*lnd;int JDS[3];xj,yj,xo,yo;float B[3][6],b1[3][6];D[3][3],E,u,t;S[3][6],s1[6][3];K[6][6];Kk[2*jds][2*jds]={{0}};float SK[2*jds][(d+1)*2]={{0}};float x1,x2,x3,y1,y2,y3,ae,*xy;FILE *fp1,*fp2,*fp3,*fp4,*fp5,*fp7;fp1=fopen("节点坐标.txt","r"); fp2=fopen("单元节点编号.txt","r");fp3=fopen("材料参数.txt","r");fp4=fopen("输出结果.txt","a"); fp5=fopen("节点总数和单元总数.txt","r");fp7=fopen("SK矩阵.txt","w"); fscanf(fp5,"%d",&nj);fscanf(fp5,"%d",&ne);xy=(float*)malloc(nj*2*sizeof(float));lnd=(int*)malloc(ne*3*sizeof(int));for(i=0;i<nj;i++)for(j=0;j<2;j++)fscanf(fp1,"%f",&xy[i*2+j]);for (i=0;i<ne;i++)for (j=0;j<3;j++)fscanf(fp2,"%d",&lnd[i*3+j]);for(n=0;n<ne;n++){printf("请输入所求的单元号:");scanf("%d",&a);i=lnd[(a-1)*3+0],j=lnd[(a-1)*3+1],k=lnd[(a-1)*3+2];JDS[0]=i;JDS[1]=j;JDS[2]=k;x1=xy[(i-1)*2+0];y1=xy[(i-1)*2+1];x2=xy[(j-1)*2+0];y2=xy[(j-1)*2+1];x3=xy[(k-1)*2+0];y3=xy[(k-1)*2+1]; ae=0.5*((x2-x1)*(y3-y1)-(x3-x1)*(y2-y1));printf("第%d单元的面积为：%7.4f\n",a,ae);fprintf(fp4,"\n第%d单元面积为：%7.4f ",a,ae);b1[0][0]=y2-y3,b1[0][2]=y3-y1;b1[0][4]=y1-y2,b1[1][1]=x3-x2;b1[1][3]=x1-x3,b1[1][5]=x2-x1;b1[2][0]=b1[1][1],b1[2][1]=b1[0][0];b1[2][2]=b1[1][3],b1[2][3]=b1[0][2];b1[2][4]=b1[1][5],b1[2][5]=b1[0][4];b1[1][0]=0,b1[0][1]=0, b1[1][2]=0,b1[0][3]=0,b1[1][4]=0,b1[0][5]=0; for(i=0;i<3;i++)for(j=0;j<6;j++)B[i][j]=(0.5/ae)*b1[i][j];fscanf(fp3,"%f %f %f",&E,&u,&t);printf("弹性模量、泊松比和厚度分别为:%7.4f %7.4f %7.4f\n",E,u,t);fprintf(fp4,"\n弹性模量、泊松比和厚度分别为:%7.4f %7.4f %7.4f",E,u,t);D[0][0]=E/(1-u*u), D[0][1]=E*u/(1-u*u), D[0][2]=0;D[1][0]=E*u/(1-u*u),D[1][1]=E/(1-u*u), D[1][2]=0;D[2][0]=0, D[2][1]=0, D[2][2]=E*(1-u)/(2*(1-u*u));for(i=0;i<3;i++)for(j=0;j<6;j++)for(k=0,S[i][j]=0;k<3;k++)S[i][j]=S[i][j]+D[i][k]*B[k][j];printf("第%d单元的应力矩阵为:\n",a);fprintf(fp4,"\n单元%d单元的应力矩阵\n",a);for(i=0;i<3;i++){for(j=0;j<6;j++){printf("%7.4f ",S[i][j]);fprintf(fp4,"%7.4f ",S[i][j]);}fprintf(fp4,"\n");printf("\n");}for(i=0;i<3;i++)for(j=0;j<6;j++)s1[j][i]=S[i][j];for(i=0;i<6;i++){for(j=0;j<6;j++)for(k=0, K[i][j]=0;k<3;k++)K[i][j]=K[i][j]+s1[i][k]*B[k][j]*ae*t;}printf("第%d单元的单元刚度阵为:\n",a);fprintf(fp4,"\n第%d单元的单元刚度阵为\n",a);for(i=0;i<6;i++){for(j=0;j<6;j++){ printf("%7.4f ",K[i][j]);fprintf(fp4,"%7.4f ",K[i][j]);}fprintf(fp4,"\n");printf("\n");}for(i=0;i<6;i++){k=i/2;xj=JDS[k]*2-1;xo=JDS[k]*2;if(i%2==0)for(j=0;j<6;j++){k=j/2;yj=JDS[k]*2-1;yo=JDS[k]*2;if(j%2==0)if(yj>=xj&&(yj-xj)<(d+1)*2)SK[xj-1][yj-xj]=SK[xj-1][yj-xj]+K[i][j];if(j%2!=0)if(yo>=xj&&(yo-xj)<(d+1)*2)SK[xj-1][yo-xj]=SK[xj-1][yo-xj]+K[i][j];}elsefor(j=0;j<6;j++){k=j/2;yj=JDS[k]*2-1;yo=JDS[k]*2;if(j%2==0)if(yj>=xo&&(yj-xo)<(d+1)*2)SK[xo-1][yj-xo]=SK[xo-1][yj-xo]+K[i][j];if(j%2!=0)if(yo>=xo&&(yo-xo)<(d+1)*2)SK[xo-1][yo-xo]=SK[xo-1][yo-xo]+K[i][j];}}}for(i=0;i<2*jds;i++){for(j=0;j<2*(d+1);j++){fprintf(fp7,"%7.4f ",SK[i][j]);}fprintf(fp7,"\n");}fclose(fp1);fclose(fp2);fclose(fp3);fclose(fp4);fclose(fp5);fclose(fp7);//fclose(fp6); }。