有限单元法作业非线性分析 程序

几何非线性大作业荷载增量法

和弧长法程序设计

一、几何非线性大作业（Newton-Raphson法）

用荷载增量法（Newton-Raphson法）编写几何非线性程序：

（1）用平面梁单元，可分析平面杆系

（2）算例：悬臂端作用弯矩。悬臂梁最终变形形成周长为悬臂梁长度的圆。

1.1 Newton-Raphson算法基本思想

图1.1 Newton-Raphson算法基本思想

1.2 悬臂梁参数

基本参数：L=2m, D=0.03m, A=7.069E-4m2, I=3.976E-08m4 ,E=2.0E11N/m2

图1.2 悬臂梁单元信息

将悬臂梁分成10个单元，如图1.2所示

2.1 MATLAB输入信息

材料信息单元信息

约束信息（0为约束，1为放松）荷载信息（FX,FY,M）

节点信息

2.2 求解过程

梁弯成圆形：理论弯矩M=EIY"=24981.944N.m ,直径为0.642m 运用ABAQUS和MATLAB进行求解对比：

图1.3 加载图

图1.4 ABAQUS变形图

图1.5 MATLAB变形曲线

ABAQUS和MATLAB变形对比，最终在理论荷载作用下都弯成了一个圆，其直径为0.64716m，与理论值相对比值为：（0.64716-0.642）/0.642=0.00804.非常接近。

2.3 加载点荷载位移曲线

图1.5 加载点Y方向的荷载位移曲线

加载点的最大竖向位移分别为1.4525m和1.45246m，相对比值（1.4525-1.45246）/1.45246=2.75395E-05。完全相同，说明MATLAB的计算结果很好。

二、几何非线性大作业（弧长法）

用弧长法编写几何非线性程序，分析荷载位移全过程曲线：

1) 用平面梁单元，可分析平面杆系结构

2) 算例

(1)受集中荷载的拱：考察拱的矢跨比、荷载位置对荷载位移曲线的影响。

(2)其他有复杂平衡路径的结构

3) 将结果与相关文献进行对比

1.1 弧长法基本思想

图2.1 弧长法基本思想

1.2 拱基本参数

拱参数：L=100m, A=0.32m2,I=1m4 ,E=1.0e7N/m2，F=-5000N，拱曲线y=5×sin(3.1415926*x/L)

将拱结构分成25个单元，如图2所示

图2.2拱单元信息

2.1 MATLAB输入信息

材料信息单元信息

约束信息（0为约束，1为放松）荷载信息（FX,FY,M）

节点信息

2.2运用ANSYS和MATLAB进行求解对比（两端铰接）ANSYS中模型:

图2.3 ANSYS模型

图2.4 MA TLAB和ANSYS变形图

2.3 加载点荷载位移曲线

图2.5 加载点荷载位移曲线

ANSYS求得的极限承载力3042.53,对应位移3.00142

MATLAB求得的极限承载力3043.8, 对应位移3.0768

相对误差分别为0.0417%,2.45%，模拟效果比较好。

2.4 拱的矢跨比a对拱荷载位移曲线的影响

不同矢跨比（1/20，3/40，1/10，3/20）下加载点的荷载位移曲线

1）MATLAB中计算拱的矢跨比a对拱荷载位移曲线的影响

图2.6 荷载位移曲线

图2.7荷载位移曲线

表1 各矢跨比下拱结构的极限荷载

参数

极值点F(N) 位移（m）最低点F(N) 位移（m）矢高

5mm 3043.8 3.0768 1765.2 7.0816 7.5mm 7623.3 4.0335 -595.82 11.21

10mm 14974 5.4026 -6408.1 14.886 20mm 39791 9.4831 -63049 30.513

从表中可以初步得出：在一定随着矢跨比的增加，拱仍然呈现跳跃失稳的形式，拱结构的极限承载能力有大幅度的提高；在最低处的承载力呈现出反向，相当于有一个拉力在阻止拱结构发生跳跃失稳，矢跨比越大，拱越不容易发生跳跃失稳。当拱的矢跨比超过一定范围后，拱将发生复杂的不同于跳跃失稳的失稳形式。

2）MATLAB与ANSYS计算结果对比

图2.8 ANSYS和MA TLAB对比荷载位移曲线

表2 各矢跨比下拱结构的极限荷载对比

参数

F(N)MAT 位移（m）F(N)ANA 位移（m）误差（%）误差（%）矢高

5mm 3043.8 3.0768 3042.53 3.00142 0.04 2.45 7.5mm 7623.3 4.0335 7624.91 3.96303 -0.02 1.75

10mm 14974 5.4026 14974.3 5.3157 0.00 1.61 20mm 39791 9.4831 39695.7 9.59955 0.24 -1.23

从图中可以看出：矢跨比在一定范围内，MATLAB与ANSYS计算的荷载位移曲线非常吻合，验证了MATLAB程序的可行性。当矢跨比为0.15时，ANSYS 中将跟踪不到失稳后复杂的平衡路径。

从表中可以得出：MATLAB与ANSYS计算中拱的极限荷载和极限荷载时所对应的位移非常接近，加载点均为顶点26。具体为：矢高5mm，荷载误差为0.04，位移误差为2.45；矢高7.5mm，荷载误差为-0.02，位移误差为1.75；矢高10mm，荷载误差为0，位移误差为-1.61；矢高20mm，荷载误差为0.24，位移误差为-1.23。实际误差相差很小，在工程允许的范围内是可以接受的。

2.5荷载位置对拱荷载位移曲线的影响

图2.9 ANSYS模型简图

1）MATLAB中计算荷载位置对拱荷载位移曲线的影响

图2.10 各加载点荷载位移曲线

表3改变加载点拱结构的极限荷载

参数

极值点F(N) 位移（m）最低点F(N) 位移（m）加载点

26 3043.8 3.0768 1765.2 7.0816

16 3570 3.1891 2258.8 6.116

11 4728 2.88 3220.5 4.7959

4 14317 1.2826 10569 1.7829

误差=100*（MATLAB-ANSYS）/ANSYS 从表中可以初步得出：随着加载点位置越靠近支座处，拱结构的极限承载能力有大幅度的提高；在最低处的承载力也大幅度提高。当加载点位置靠近支座

时，拱的承载力增加幅度最大，拱的稳定性很强，不容易发生失稳。

2）MATLAB与ANSYS计算结果对比

图2.11 ANSYS和MATLAB对比荷载位移曲线

表4各加载点拱结构的极限荷载

参数

F(N)MAT 位移（m）F(N)ANA 位移（m）误差（%）误差（%）矢高

26 3043.8 3.0768 3042.53 3.00142 0.04 2.45

16 3570 3.1891 3569.73 3.24865 0.01 -1.87

11 4728 2.88 4728.71 2.91862 -0.02 -1.34

4 14317 1.2826 14324.8 1.29764 -0.0

5 -1.17

误差=100*（MATLAB-ANSYS）/ANSYS 从图中可以看出：MATLAB与ANSYS计算的荷载位移曲线非常吻合，验证了MATLAB程序的可行性。从表中可以得出：MATLAB与ANSYS计算中拱的极限荷载和极限荷载时所对应的位移非常接近。具体为：26加载点，荷载误差为0.04，位移误差为2.45；16加载点，荷载误差为0.01，位移误差为-1.87；11加载点，荷载误差为-0.02，位移误差为-1.34；4加载点，荷载误差为-0.05，位移误差为-1.17。实际误差相差很小，在工程允许的范围内是可以接受的。

2.6 两端铰接和固接对拱荷载位移曲线的影响

矢高为5mm 时，拱两端为固接和铰接时的荷载位移曲线如下：

图2.12 ANSYS 和MA TLAB 固接和铰接的荷载位移曲线

从图中可以看出：拱的边界条件对其的失稳形式有很大影响。两端固接拱的稳定性明显优于两端铰接拱，承载能力也大幅度提高。固接拱不会发生跳跃失稳的形式，刚度在初始阶段会减小，待到达一定程度后刚度又会增加。而两端铰接拱会发生跳跃失稳的形式。

2.7参数m 对拱失稳形式的影响

文献中给出:m 是一个由拱截面在竖向平面内的回转半径r 和拱的初始矢高h 无确定的无量纲参数。

2

242()I

r m h Ah

==

当m>=1 时，扁拱不会出现跳跃屈曲, 当0

2.13 不同m值加载点的荷载位移曲线

2.14 不同m值变形后拱曲线

从MATLAB的计算结果中可以验证：不同的m系数对应拱不同的失稳形式。

当m>=1 时，扁拱不会出现跳跃屈曲，当0

2.8 具有复杂失稳形式的拱

铰支圆拱

该结构及其几何参数、物理性质均示于图4a 中。中心受集中荷载。这个结构是研究分歧问题的经典题目。将半跨结构划分为8个单元, 得到图4b的基本路径和分歧路径, 并与JChreseielewski和Rsehmiot的结果进行了比较。本文对此结构也进行了缺陷分析。

拱的基本参数：L=254cm，R=381cm，I=41.62cm4，A=6.45cm2，E=6898kN/cm2。

文献中的计算结果。

采用MATLAB和ANSYS对其进行求解，得到其荷载位移曲线：

图2.15 ABAQUS模型

图2.16 ABAQUS变形图

图2.17 ANSYS、MA TLAB、ABAQUS加载点荷载位移曲线从MATLAB、ANSYS、ABAQUS计算的荷载位移曲线可以看出：两者的荷载位移曲线基本吻合。MATLAB的计算结果可以看出在后期其承载能力会有较大提高，与文献中的荷载位移曲线趋势相同，所以验证出程序的可靠性。ABAQUS不能模拟出后续段曲线也许是单元划分过少。

图2.18 MATLAB加载点荷载位移曲线

第二次极值点会超过第一次极值点所对应的荷载，与文献一致，荷载点也接近。

加入初始缺陷：L/1000, L/2000初始缺陷后观察加载点的荷载位移曲线变化趋

势。

图2.19 ANSYS加入初始缺陷后的加载点荷载位移曲线

2.20 初始缺陷为0.0001时的荷载位移曲线

加入初始缺陷后，拱的极限承载能力明显降低。且失稳形式也发生了变化，初始缺陷的大小对其荷载位移曲线有明显影响。所以在工程设计中应考虑结构或构件的初始缺陷，使结构的设计更加合理，安全。

为了研究初始缺陷对拱失稳路径的影响，应用ABAQUS和ANSYS以及MATLAB中加水平力模拟拱结构初始缺陷下的荷载位移曲线。

为了探究ABAQUS和ANSYS计算结果，取其前三阶模态进行对比分析：

2.21一阶屈曲模态

ABAQUS和MATLAB中的一阶屈曲系数为0.55884和0.564512，对应的屈曲荷载为74325.72N 和75080.096N。

2.22 二阶屈曲模态

ABAQUS和MA TLAB中的二阶屈曲系数为1.2259和1.253，对应的屈曲163044.7N 和

非线性大作业

非线性大作业 -标准化文件发布号：（9456-EUATWK-MWUB-WUNN-INNUL-DDQTY-KII

工程结构非线性分析 学院： 姓名： 学号： 指导教师：

目录 1、预应力混凝土梁截面非线性 0 1.1 材料的本构关系 0 1.2 平截面假定 (2) 1.3 预应力筋作用下截面初应变的求解 (3) 2、预应力混凝土梁构件的非线性 (4) 2.1 构件弯曲的一般理论 (4) 2.2 共轭梁分析法 (4) 2.3 预应力钢筋混凝土梁非线性分析的数值法 (5) 3、算例分析 (7) 3.1 试验梁简介 (7) 3.2 截面非线性与构件非线性分析程序编制 (8) 3.3 试验结果验证 (9) 3.4 结果分析 (11) 参考文献 (12) 附录 (13)

作业2：预应力混凝土梁的非线性全过程分析 要求： 1.阐述预应力混凝土梁截面和构件非线性全过程分析的理论背景； 2.编制相应的截面和构件非线性分析程序，给出具体算例分析结果，方法及程序的适用性必须有试验结果的验证。

1、预应力混凝土梁截面非线性 1.1 材料的本构关系 1.1.1 混凝土本构关系 混凝土受压采用Rush 建议的应力—应变曲线，如图1-1所示。 0cu f c σ 图1-1 混凝土受压应力-应变曲线 000[1(1)]n c c c c c c cu f f εεεεεεε ? -- 0≤≤?σ=? ? <≤? 式中c σ——对应于混凝土应变为c ε时的混凝土压应力； c f ——混凝土抗压强度标准值； cu ε——正截面处于非均匀受压时的混凝土极限压应变， 50.0033(50)10cu cu f ε-=--?，当0.0033cu ε>时，取为0.0033； 0ε——受压峰值应变，500.0020.5(50)10cu f ε-=+-?，当00.002ε<时，取为0.002； n ——系数，,1 2(50)60 cu k n f =- -，当 2.0n >时，取为2.0。 为计算方便，混凝土受拉应力-应变曲线采用线性式，如图1-2所示。

非线性分析作业讲义

学院：材料科学与工程学院专业：材料工程 姓名：飞学号：1125 作业： 找出几个所在专业研究领域的重要而且有研究价值的非线性问题及其模型，要求写出相应的模型方程及其所涉及的变量参数涵义，并列举出研究该模型的主要研究现状。（不少于3种） 举例1：材料力学领域的非线性问题 非线性本构和非线性本构复合材料 1.1 研究非线性本构模型的意义 从力学的角度来看，C/SiC复合材料属于准脆性的各向异性材料。以碳纤维、热解碳界面和SiC基体三种典型组分构成的C/SiC复合材料为例，相对于脆性的单质陶瓷，该材料具有较好的韧性。主要原因是在机械载荷作用下，材料内部存在如前所述的基体开裂、界面脱粘和滑移、纤维断裂和拔出等多种能量耗散机制。虽然这些细观损伤模式有别于金属的屈服机理，但是材料表现出类似的弹塑性-损伤力学行为。图1-1为C/SiC复合材料在沿轴向拉伸加卸载条件下的典型应力-应变曲线，从图中可看出：材料的线弹性极限较低，通常为20MPa左右；当应力水平超过弹性极限之后，材料的弹性模量（E0）开始减小，同时产生类似于不可回复的残余应变，卸载-重加载过程中应力-应变曲线形成迟滞环，且迟滞环的宽度随卸载点应力的增大而不断增大。该材料的剪切应力-应变关系也有类似的特征。由此易知，在对C/SiC复合材料的应力-应变关系进行分析描述时，传统的线弹性本构模型已经不再胜任；而如果仅在线弹性范围内使用该材料，则不能充分发挥出材料的力学性能，安全裕度过大，与航空航天器追求减重的目标不符。因此需要充分了解该材料的非线性力学行为，特别是其内部的损伤机理与特性，并为其建立合适的非线性本构模型。

图1-1 C/SiC复合材料的典型拉伸加/卸载应力-应变曲线 建立非线性本构模型的一个重要作用是辅助C/SiC复合材料的结构优化设计。如前所述，目前C/SiC复合材料已经开始逐步在航空航天器结构上使用，轻质、可重复使用等特性有助于提高飞行器的性能，并降低寿命周期内的使用和维护成本，但是这类材料仍然存在造价高的缺点。例如，德国DLR为X-38 V201飞行器提供的全C/SiC复合材料襟翼的尺寸约为1.4m×1.6m，重68公斤，造价高达2千万美元。这是由材料制备工艺的特点决定的。以较为成熟的等温CVI 工艺为例，该工艺具有能够制备出高纯度的基体、可用于一定厚度构件的近尺寸成型等诸多优点，但是为防止沉积的基体太快地封堵预制体孔隙通道，需要在相对缓慢的沉积速率下进行，因此材料的制备周期长，通常需要几周或数百小时的时间，而且化学反应过程中生成的HCl等副产物对设备有腐蚀作用，导致制备成本偏高，限制了材料的推广应用。因此，为C/SiC复合材料建立合适的本构模型，在结构设计阶段将本构模型与商业有限元软件结合，准确计算和结构在不同受载条件下的应力状态并预测其承载能力，有助于结构的优化设计，同时省去或减少大量的试件制备和测试过程，从而降低热结构的研发成本。国内已经对C/SiC 的损伤机理和本构模型开展了一些研究工作。潘文革等人对二维和三维编织C/SiC复合材料在单轴拉伸载荷下的损伤演化进行了试验研究，通过分析声发射事件数和相对能量等参数，发现两种材料的拉伸损伤过程大致分为初始损伤阶段、过渡阶段、损伤加速和快速断裂阶段；杨成鹏等人对二维编织C/SiC复合材料单轴拉伸非线性力学行为进行了试验研究，通过循环加卸载试验方法，获得了材料的残余应变和卸载模量随拉伸应力的变化关系，并建立了基于剪滞理论的细观损伤力学模型；陶永强等人将二维编织结构简化成正交铺层和纤维束波动部分的组合，采用了Curtin和Ahn提出的基体随机开裂、纤维随机断裂的统计分布理论以及体积平均方法，预测了二维编织C/SiC复合材料的应力-应变关系。此

有限元基础知识归纳

有限元知识点归纳 1.、有限元解的特点、原因？ 答：有限元解一般偏小，即位移解下限性 原因：单元原是连续体的一部分，具有无限多个自由度。在假定了单元的位移函数后，自由度限制为只有以节点位移表示的有限自由度，即位移函数对单元的变形进行了约束和限制，使单元的刚度较实际连续体加强了，因此，连续体的整体刚度随之增加，离散后的刚度较实际的刚度K为大，因此求得的位移近似解总体上将小于精确解。 2、形函数收敛准则（写出某种单元的形函数，并讨论收敛性）P49 (1)在节点i处N i=1,其它节点N i=0; (2)在单元之间，必须使由其定义的未知量连续； (3)应包含完全一次多项式； (4)应满足∑Ni=1 以上条件是使单元满足收敛条件所必须得。可以推证，由满足以上条件的形函数所建单元是完备协调的单元，所以一定是收敛的。 4、等参元的概念、特点、用时注意什么？（王勖成P131） 答：等参元—为了将局部坐标中几何形状规则的单元转换成总体（笛卡尔）坐标中的几何形状扭曲的单元，以满足对一般形状求解域进行离散化的需要，必须建立一个坐标变换。即： 为建立上述的变换，最方便的方法是将上式表示成插值函数的形式，即： 其中m是用以进行坐标变换的单元节点数，xi,yi,zi是这些结点在总体（笛卡尔）坐标内的坐标值，Ni’称为形状函数，实际上它也是局部坐标表示的插值函数。称前者为母单元，后者为子单元。 还可以看到坐标变换关系式和函数插值表示式：在形式上是相同的。如果坐标变换和函数插值采用相同的结点，并且采用相同的插值函数，即m=n，Ni’=Ni,则称这种变换为等参变换。 5、单元离散？P42 答：离散化既是将连续体用假想的线或面分割成有限个部分，各部分之间用有限个点相连。每个部分称为一个单元，连接点称为结点。对于平面问题，最简单、最常用的离散方式是将其分解成有限个三角形单元，单元之间在三角形顶点上相连。这种单元称为常应变三角形单元。常用的单元离散有三节点三角形单元、六节点三角形单元、四节点四边形单元、八节点四边形单元以及等参元。 6、数值积分，阶次选择的基本要求？ 答：通常是选用高斯积分 积分阶次的选择—采用数值积分代替精确积分时，积分阶数的选取应适当，因为它直接影响计算精度，计算工作量。选择时主要从两方面考虑。一是要保证积分的精度，不损失收敛性；二是要避免引起结构总刚度矩阵的奇异性，导致计算的失败。

《有限单元法》编程作业

湖南大学 《有限单元法》编程大作业 专业：土木工程 姓名： 学号： 2013年12月

目录 程序作业题目： (3) 1、程序编制总说明 (3) 2、Matlab程序编制流程图 (3) 3、程序主要标示符及变量说明 (4) 4、理论基础和求解过程 (5) 4.1、构造插值函数 (5) 4.2位移插值函数及应变应力求解 (5) 5.程序的验证 (6) 附录：程序代码 (15)

程序作业题目： 完成一个包含以下所列部分的完整的有限元程序( Project) 须提供如下内容的文字材料（1500字以上）： ①程序编制说明； ②方法的基本理论和基本公式； ③程序功能说明； ④程序所用主要标识符说明及主要流程框图； ⑤ 1～3 个考题：考题来源、输出结果、与他人成果的对比结果（误差百分比）； ⑥对程序的评价和结论（包括正确性、适用范围、优缺点及其他心得等）。 须提供源程序、可执行程序和算例的电子文档或文字材料。选题可根据各自的论文选题等决定。 1、程序编制总说明 a.该程序采用平面三角形等参单元，能解决弹性力学的平面应力、平面应变问题。 b.能计算单元受集中力的作用。 c.能计算结点的位移和单元应力。 d.考题计算结果与理论计算结果比较，并给出误差分析。 e.程序采用MATLAB R2008a编制而成。 2、Matlab程序编制流程图

图1 整个程序流程图 3、程序主要标示符及变量说明 1、变量说明： Node ------- 节点定义 gElement ---- 单元定义 gMaterial --- 材料定义，包括弹性模量，泊松比和厚度 gBC1 -------- 约束条件 gNF --------- 集中力 gk------------总刚 gDelta-------结点位移 输入结构控制参数 输入其它数据 形成整体刚度阵 引入支承条件 解方程，输出位移 求应力，输出应力 形成节点荷载向量 开始 结束 1 单元面积 求弹性矩阵 单元刚度矩阵 位移-应变矩阵 6 7 8 9 10 2 3 4 5

实例matlab-非线性规划-作业

实例matlab-非线性规划-作业

现代设计方法-工程优化理论、方法与设计 姓名 学号 班级 研 问题 ： 某厂向用户提供发动机，合同规定，第一、二、三季度末分别交货40台、60台、80台。每季度的生产费用为 （元），其中x 是该季生产的台数。若交货后有剩余，可用于下季度交货，但需支付存储费，每台每季度c 元。已知工厂每季度最大生产能力为100台，第一季度开始时无存货，设a=50、b=0.2、c=4，问工厂应如何安排生产计划，才能既满足合同又使总费用最低。讨论a 、b 、c 变化对计划的影响，并作出合理的解释。 问题的分析和假设： 问题分析：本题是一个有约束条件的二次规划问题。决策变量是工厂每季度生产的台数，目标函数是总费用（包括生产费用和存储费）。约束条件是生产合同，生产能力的限制。在这些条件下需要如何安排生产计划，才能既满足合同又使总费用最低。 问题假设： 1、工厂最大生产能力不会发生变化； 2、合同不会发生变更； 3、第一季度开始时工厂无存货； 4、生产总量达到180台时，不在进行生产； 5、工厂生产处的发动机质量有保证，不考虑退货等因素； 6、不考虑产品运输费用是否有厂家承担等和生产无关的因素。 符号规定： x1——第一季度生产的台数； x2——第二季度生产的台数； 180-x1-x2——第三季度生产的台数； y1——第一季度总费用； y2——第二季度总费用； y3——第三季度总费用； y ——总费用（包括生产费用和存储费）。 ()2bx ax x f +=

建模： 1、第一、二、三季度末分别交货40台、60台、80台； 2、每季度的生产费用为 （元）； 3、每季度生产数量满足40 ≤x1≤100，0≤x2≤100，100≤x1+x2 ≤180； 4、要求总费用最低，这是一个目标规划模型。 目标函数: y1 2111x b x a Z ?+?= y2()4012222-?+?+?=x c x b x a Z y3()()()10018018021221213 -+?+--?+--?=x x c x x b x x a Z y x x x x x x Z Z Z Z 68644.04.04.0149201 212221321--+++=++= 40≤x1≤100 0≤x2≤100 100≤x1+x2≤180 ()2 bx ax x f +=

有限单元法作业

ANSYS静力分析论文 学院：能源与动力工程学院 姓名：王立伟 班级：热能1003 学号：10110303

1 前言 随着现代工业的不断发展，人们对产品质量的要求逐步提高，传统的产品设计技术目前已远远不能满足产品的功能和市场的要求。而现代设计技术是以电子计算机为手段，以网络为基础，建立在现在管理之上，运用工程设计的新理论、新方法，实现计算机结果最优化，设计过程高效化的设计技术，它是传统设计技术的延伸和发展，它使传统设计技术发生了质的飞跃。 有限元法已成为非常普及的数字化分析方法，国际上已发布了众多的有限元分析软件，因此，甚至可以说只要你能够进行工程设计和画图，就可以进行有限元分析。下面对实际工程问题简单的介绍一下机械优化设计方法的应用。 2 有限元设计部分 2.1 问题阐述 如图2-1所示为对称涵洞模型的一半【1】，将其作为平面应变问题，其顶部的路面作用有均布载荷6000N/m2。基于有限元分析软件，对该问题进行力学分析。其中定义模型的弹性模量E=210GPa，u=0.3 图2-1

2.2 近似与假设 本题的分析为平面应变问题。压力载荷只作用在X-Y平面上。近似操作是使用固体模型来构造2-D模型并利用节点和单元将其自动划分网格。 2.3 主要分析思想 为了能对有限元分析有更好的掌握，本次课设应用ANSYS对此问题进行交互式求解。求解过程分别采用PLANE42、PLANE183、PLANE2单元类型进行网格的划分，PLANE2类型的单元分别用0.2和0.05的单元边长进行网格划分。并对各种划分单元的结果和分析后的结果进行比较。 2.4 ANSYS的求解过程 2.4.1 进入ANSYS 在D盘建立一文件夹，文件名为ansys。然后运行程序→Ansys 10.0 →Configure ANSYS Products →file Management →select Working Directory: D:\ansys，input job name:Handong→Run 2.4.2 设置不显示日期和时间 Utility Menu→PlotCtrls→Window Controls→Window Options→DATE DATE/TIME display：NO DATE or TIME 2.4.3 设置计算类型 ANSYS Main Menu→Preferences…→Select Structural→OK 2.4.4 设置单元类型 Main Menu→Preprocessor→Element Type→Add/Edit/Delete→Add→Structural Solid→Quad 4-Noded 42（PLANE82）→OK→Options→在第一个（K3后的）下拉列表中选择Plane stress项→OK→Close。 2.4.5 定义材料参数 ANSYS Main Menu →Preprocessor→Material Props→Material Models→Material Models Available→Structural (双击打开子菜单) →Linear(双击) →Elastic(双击)>Isotropic(双击) →EX:210e9,PRXY: 0.3→OK 关闭材料定义菜单 2.4.6 生成几何模型 Step1生成4个关键点： ANSYS Main Menu →Preprocessor→Modeling→Creat→Keypoints →In Active CS→按次序输入3个特征点，方式为：只在X,Y,Z的3个空格内填入点的坐标，每完成一个点的输入，用Apply结束，3个特征点坐标为1(0,0,0)，2（6，0,0），3(3.2，5.5,0)和4（0,5.5,0）→OK Step2建立4条线： ANSYS Main Menu→Preprocessor →Modeling→Create→Lines→Lines→In Active Coord→用鼠标选择关键点1和2形成L1，选择关键点2和3形成L2，选择关键点3和4形成L3→选择关键点4和1形成L4→OK Step3 创建面： ANSYS Main Menu→Preprocessor→Modeling→Create→Areas→Arbitrary →By L ines→鼠标单击选择4条线→OK

数值分析大作业三 四 五 六 七

大作业 三 1. 给定初值 0x 及容许误差 ，编制牛顿法解方程f (x )=0的通用程序. 解：Matlab 程序如下： 函数m 文件：fu.m function Fu=fu(x) Fu=x^3/3-x; end 函数m 文件：dfu.m function Fu=dfu(x) Fu=x^2-1; end 用Newton 法求根的通用程序Newton.m clear; x0=input('请输入初值x0:'); ep=input('请输入容许误差:');

flag=1; while flag==1 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0)

while flag==1 sigma=k*eps; x0=sigma; k=k+1; m=0; flag1=1; while flag1==1 && m<=10^3 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0)=ep flag=0;

end end fprintf('最大的sigma 值为：%f\n',sigma); 2．求下列方程的非零根 5130.6651()ln 05130.665114000.0918 x x f x x +?? =-= ?-???解： Matlab 程序为： （1）主程序 clear clc format long x0=765; N=100; errorlim=10^(-5); x=x0-f(x0)/subs(df(),x0); n=1;

结构非线性作业参考-方志老师

《工程结构非线性》作业 学院：土木工程学院 专业：桥梁工程 姓名：刘万事 学号：S140110021 教师：方志（教授）

结构非线性作业一 （1） 求出荷载—柱中点侧移的解析解及第一类失稳荷载； （2） 以具体的实例给出几何非线性效应的数值解（可用有限元程序计算,但应给出有限元程 序理论背景的详细描述），并与解析解结果对比； （3） 对结构几何非线性和稳定的关系进行讨论。 1、偏压柱的跨中最大挠度的解析解 图1 计算简图 1.1跨中弯矩为： ()M P e y =+ （1） 1.2由材料力学中梁挠曲线的近似微分方程可以得到： 22d y M dx EI =- 将（1）式代入其中得 ''()P e y P Pe y y EI EI EI +=-=+ 解微分方程得： []csc sin ()csc sin sin ()sin csc l x x y e l l x e l x l l e l x x l e ααααααα-??? ?=--+-????????=-+-

其中α= 1.3 求跨中侧移：当2l x = 时 max 2sin csc (sec 1)22 l l y e l e e ααα=-=- 2、用有限元软件ansys 建立题中所给的弯压柱的力学模型，并计算跨中最大挠度 2.1 给出一个实例： 假设题中所给弯压柱所受荷载P=10KN, 偏心距e=0.05m ，柱高为L=5m ，采用屈服 强度为345MP 的钢材，弹性模量E=2.06x105 MP, 柱的截面尺寸如所示： 图1 计算截面 2.2 确定材料的本构模型 采用韩林海（2007）中的二次塑性流模型来模拟钢材， 其应力-应变关系曲线,分为弹性段(Oa)、弹塑性段(ab)、塑性段 (bc)、强化段(cd)和二次塑流(de)等五个阶段，如图1所示。图1中的点划线为钢材实际的应力-应变关系曲线，实线所示为简化的应力-应变关系曲线，模型的数学表达式如式(3-1)。其中： e1e3e1e2e e1s y e 100,10,5.1,/8.0εεεεεεε====E f ； f p 、f y 和 f u 分别为钢材的比例极限、屈服极限和抗拉强度极限。

电子科技大学非线性系统作业

非线性系统理论分析及其应用 XXX （1．电子科技大学，XXXXX学院，XXXXXXX） Theoretical Analysis Of nonlinear Systems And Its Applications XXXXXXXXXX (University of Electronic Science and Technology of China，School of Energy Science and Engineering,XXXXXXXXXXXXXXXX) 摘要：本文通过通过对非线性系统的原理，分类，性质等做了细致的分析，并重点介绍了非线性系统在电力系统，自行车自动控制等方面的应用，得出非线性系统在控制领域的重要地位。 关键词：非线性；原理；应用 ABSTRACT: In this paper, through the principle of non-linear systems, classification, properties, and so do a detailed analysis and focuses on the application of nonlinear systems in the power system, automatic control and other aspects of the bike, draw an important role in the control field of nonlinear systems . KEY WORDS：Nonlinear; principle; application 1 非线性系统的原理 非线性系统是状态变量和输出变量对于所有可能的输入变量和初始状态都满足叠加原理的系统。一个由线性元部件所组成的系统必是线性系统。但是，相反的命题在某些情况下可能不成立。线性系统的状态变量(或输出变量)与输入变量间的因果关系可用一组线性微分方程或差分方程来描述，这种方程称为系统的数学模型。 1.1 非线性与线性概述 线性，指量与量之间按比例、成直线的关系，在空间和时间上代表规则和光滑的运动；而非线性则指不按比例、不成直线的关系，代表不规则的运动和突变。如问：两个眼睛的视敏度是一个眼睛的几倍？很容易想到的是两倍，可实际是6－10倍！这就是非线性：1+1不等于2。激光的生成就是非线性的！当外加电压较小时，激光器犹如普通电灯，光向四面八方散射；而当外加电压达到某一定值时，会突然出现一种全新现象：受激原子好像听到“向右看齐”的命令，发射出相位和方向都一致的单色光，就是激光。非线性的特点是：横断各个专业，渗透各个领域，几乎可以说是：“无处不在时时有。”如：天体运动存在混沌；电、光与声波的振荡，会突陷混沌；地磁场在400万年间，方向突变16次，也是由于混沌。甚至人类自己，原来都是非线性的：与传统的想法相反，健康人的脑电图和心脏跳动并不是规则的，而是混沌的，混沌正是生命力的表现，混沌系统对外界的刺激反应，比非混沌系统快。由此可见，非线性就在我们身边，躲也躲不掉了。 1.2 非线性与线性的比较 定性地说，线性关系只有一种，而非线性关系则千变万化，不胜枚举。线性是非线性的特例，它是简单的比例关系，各部分的贡献是相互独立的；而非线性是对这种简单关系的偏离，各部分之间彼此影响，发生耦合作用，这是产生非线性问题的复杂性和多样性的根本原因。正因为如此，非线性系统中各种因素的独立性就丧失了：整体不等于部分之和，叠加原理失效，非线性方程的两个解之和不再是原方程的解。因此，对于非线性问题只能具体问题具体分析。 线性与非线性现象的区别一般还有以下特征：（1）在运动形式上，线性现象一般表现为时空中的平滑运动，并可用性能良好的函数关系表示，而非线性现象则表现为从规则运动向不规则运动的转化和跃变； （2）线性系统对外界影响的响应平缓、光滑，而非线性系统中参数的极微小变动，在一些关节点上，可以引起系统运动形式的定性改变。在自然界和人类社会中大量存在的相互作用都是非线性的，线性作用只不过是非线性作用在一定条件下的近似。 1.3 非线性系统分类

有限元法大作业

有限元法大作业 一平面刚架的程序 用Visual C++编制的平面刚架的源程序如下： ///////////////////////////////////////////////////////程序开始////////////////////////////////////////////////////////////////// #include"iostream.h" #include"math.h" #include"stdlib.h" #include"conio.h" //***************** //声明必要变量 //***************** #define PI 3.141592654 int NE; //单元数 int NJ; //节点数 int NZ; //支承数 int NPJ; //有节点载荷作用的节点数 int NPF; //非节点载荷数 int HZ; //载荷码 int E; //单元码 int fangchengshu; double F[303]; //各节点等效总载荷数值 int dym_jdm[100][2]; //单元码对应的节点码:dym_jdm[][0], dym_jdm[][1]分为前后节点总码 int zhichengweizhi[300]; //记录支持节点作用点的数组 int fjzhzuoyongdanyuan[100]; //非节点载荷作用单元 int fjzhleixing[100]; //非节点载荷类型:1-均布,2-垂直集中,3-平行集中,4-力偶,5-角度集中 double fjzhzhi[100]; //非节点载荷的值 double fjzhzuoyongdian[100]; //非节点载荷在各竿的作用点 double fjzhjiaodu[100]; //非节点载荷作用角度 int jdzhzuoyongdian[100]; //节点载荷作用的节点数组 double jiedianzaihe[101][3];//节点载荷值,其jiedianzaihe[][0]-- jiedianzaihe[][2]分别为U, V, M double zhengtigangdu[303][303]; //整体刚度数组 double changdu[100]; //各单元竿长数组 double jiaodu[100]; //各单元角度数组 double tanxingmoliang[100]; //各单元弹性模量数组 double J_moliang[100]; //各单元J模量数组 double mianji[100]; //各单元面积数组 double weiyi[303]; //记录各个节点位移的数组 double dy_weiyi[100][6]; //各个单元在局部坐标系中的位移数组dy_weiyi[i][0]-dyweiyi[i][6]分别为第i+1单元的u1,v1,@1,u2,v2,@2 double dy_neili[100][6];//各个单元在局部坐标系中的固端内力dy_weiyi[i][0]-dyweiyi[i][6]分别为第i+1单元的U1,V1,M1,U2,V2,M2 double gan_neili[100][6];//各个单元的竿端内力数组,gan_neili[i][6]表示第i+1单元的6内力. //*******************

数值计算大作业——刘剖析

课程设计 课程名称：数值分析 设计题目：数值计算大作业 学号：S315070064 姓名：刘峰 完成时间：2015年10月25日

题目一、非线性方程求根 1.题目 假设人口随时间和当时人口数目成比例连续增长，在此假设下人口在短期内的增长建立数学模型。 （1）如果令()N t 表示在t 时刻的人口数目，β表示固定的人口出生率，则人口数目满足微分方程 () ()dN t N t dt β=，此方程的解为0()=t N t N e β； （2）如果允许移民移入且速率为恒定的v ，则微分方程变成() ()dN t N t v dt β=+， 此方程的解为0()=+ (1)t t v N t N e e βββ -； 假设某地区初始有1000000人，在第一年有435000人移入，又假设在第一年年底该地区人口数量1564000人，试通过下面的方程确定人口出生率β，精确到410-；且通过这个数值来预测第二年年末的人口数，假设移民速度v 保持不变。 435000 1564000=1000000(1)e e βββ + - 2.数学原理 采用牛顿迭代法，牛顿迭代法的数学原理是，对于方程0)(=x f ，如果)(x f 是线性函数，则它的求根是很容易的，牛顿迭代法实质上是一种线性化方法，其基本思想是将非线性方程0)(=x f 逐步归结为某种线性方程来求解。 设已知方程0)(=x f 有近似根k x (假定0)(≠'x f )，将函数)(x f 在点k x 进行泰 勒展开，有 . ))(()()(???+-'+≈k k k x x x f x f x f 于是方程0)(=x f 可近似地表示为 ))(()(=-'+k k x x x f x f 这是个线性方程，记其根为1k x +，则1k x +的计算公式为 )() (1k k k k x f x f x x '- ==+，,,2,1,0???=k

大连理工大学非线性分析报告第三次作业

1）分叉图为： 程序为： muv=0:0.002:3; %%%分叉参数 m=length(muv); for k=1:m mu=muv(k); n=1000;x=zeros(n,1);x(1)=0.1; for kk=2:n x(kk)=mu*x(kk-1)*(1-x(kk-1)*x(kk-1)); %%%映射end figure(5) plot(zeros(50,1)+mu,x(301:350),'r.'); hold on; xlabel('a');ylabel('x_n');title('分叉图') grid on end

Lyapunov 图 程序为： n=400; xn=zeros(1,n); aa=2.4:0.01:4; N=1; hold on ;box on ;xlim([min(aa),max(aa)]); XL(1)=ylabel('\itx'); for a=aa; x=0.1; for q=1:80; x=a*x*(1-x); end s=0; for q=1:n; xn(q)=x; df=a-2*a*x; s=s+log(abs(df)); x=a*x*(1-x); end L(N)=s/n; N=N+1; 2.4 2.6 2.83 3.2 3.4 3.6 3.84 a

a,pause(0.01) end plot(aa,L); hold on ;box on ; grid on; xlim([min(aa),max(aa)]); 2 求解方程131n n n n n n x y y bx dy y ++=???=-+-?? 可得到固定点 ()()11,0,0x y = 当10d b -->时， ( )22,y x =，( )(33,y x =， 求得特征值表达式为 ， 利用相图，观察strange attractor ： 当b=0.2，d=2.5或2.65时，不存在奇异吸引子； 当b=0.2，d=2.77时，存在奇异吸引子，奇异吸引子图为： 234)3(222d y b d y ----±= λ

非线性分析作业第2次(硕士博士非线性分析)

1. For the following dynamical systems 1）'' 30x x x ++= 2）' '2(1),(1)3x x x xy y y y xy =--=-- a) Find all fixed points and classify them. b) Sketch the phase space portrait. Solution for 1）：'' 3 0x x x ++= Set 121 ,y x y y '==. Then, the equation becomes to ， 1 23 2 11y y y y y '=??'=--? Set vector variable z, we can write ()z f z =, where 12y z y ?? =???? 213211()y y f z f y y y ????==????--???? There is only fixed point 00z ?? =???? The Jacobian matrix 2 10 1310Df y ??=??--?? Jacobian matrix for linearized system at the fixed point, () 0110Df y ?? =?? -?? Eigenvalues for this system are 12i λ=±, so they have zero real part and the method of linearization cannot decide about the stability.

有限单元法

有限单元法 有限元方法的基础是变分原理和加权余量法，其基本求解思想是把计算域划分为有限个互不重叠的单元，在每个单元内，选择一些合适的节点作为求解函数的插值点，将微分方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式，借助于变分原理或加权余量法，将微分方程离散求解。采用不同的权函数和插值函数形式，便构成不同的有限元方法。有限元方法最早应用于结构力学，后来随着计算机的发展慢慢用于流体力学的数值模拟。在有限元方法中，把计算域离散剖分为有限个互不重叠且相互连接的单元，在每个单元内选择基函数，用单元基函数的线形组合来逼近单元中的真解，整个计算域上总体的基函数可以看为由每个单元基函数组成的，则整个计算域内的解可以看作是由所有单元上的近似解构成。在河道数值模拟中，常见的有限元计算方法是由变分法和加权余量法发展而来的里兹法和伽辽金法、最小二乘法等。根据所采用的权函数和插值函数的不同，有限元方法也分为多种计算格式。从权函数的选择来说，有配置法、矩量法、最小二乘法和伽辽金法，从计算单元网格的形状来划分，有三角形网格、四边形网格和多边形网格，从插值函数的精度来划分，又分为线性插值函数和高次插值函数等。不同的组合同样构成不同的有限元计算格式。对于权函数，伽辽金(Galerkin)法是将权函数取为逼近函数中的基函数；最小二乘法是令权函数等于余量本身，而内积的极小值则为对代求系数的平方误差最小；在配置法中，先在计算域内选取N个配置点。令近似解在选定的N个配置点上严格满足微分方程，即在配置点上令方程余量为0。插值函数一般由不同次幂的多项式组成，但也有采用三角函数或指数函数组成的乘积表示，但最常用的多项式插值函数。有限元插值函数分为两大类，一类只要求插值多项式本身在插值点取已知值，称为拉格朗日(Lagrange)多项式插值；另一种不仅要求插值多项式本身，还要求它的导数值在插值点取已知值，称为哈密特(Hermite)多项式插值。单元坐标有笛卡尔直角坐标系和无因次自然坐标，有对称和不对称等。常采用的无因次坐标是一种局部坐标系，它的定义取决于单元的几何形状，一维看作长度比，二维看作面积比，三维看作体积比。在二维有限元中，三角形单元应用的最早，近来四边形等参元的应用也越来越广。对于二维三角形和四边形电源单元，常采用的插值函数为有Lagrange插值直角坐标系中的线性插值函数及二阶或更高阶插值函数、面积坐标系中的线性插值函数、二阶或更高阶插值函数等。 对于有限元方法，其基本思路和解题步骤可归纳为 (1)建立积分方程，根据变分原理或方程余量与权函数正交化原理，建立与微分方程初边值问题等价的积分表达式，这是有限元法的出发点。 (2)区域单元剖分，根据求解区域的形状及实际问题的物理特点，将区域剖分为若干相互连接、不重叠的单元。区域单元划分是采用有限元方法的前期准备工作，这部分工作量比较大，除了给计算单元和节点进行编号和确定相互之间的关系之外，还要表示节点的位置坐标，同时还需要列出自然边界和本质边界的节点序号和相应的边界值。 (3)确定单元基函数，根据单元中节点数目及对近似解精度的要求，选择满足一定插值条件的插值函数作为单元基函数。有限元方法中的基函数是在单元中选取的，由于各单元具有规则的几何形状，在选取基函数时可遵循一定的法则。 (4)单元分析：将各个单元中的求解函数用单元基函数的线性组合表达式进行逼近；再将近似函数代入积分方程，并对单元区域进行积分，可获得含有待定系数(即单元中各节点的参数值)的代数方程组，称为单元有限元方程。 (5)总体合成：在得出单元有限元方程之后，将区域中所有单元有限元方程按一定法则进行累加，形成总体有限元方程。

非线性系统分析习题

非线性系统分析习题

第2章 2-1 电路如题图2-1所示，若11tanh 2u i =，2 23 22i i +=ψ，33ln u q =，试讨论对下列各组 变量：（1）2i 和3u ；（2）2i 和3q ；（3）2ψ和3u ；（4）2ψ和3q ；是否存在标准形式的状态方程？若存在，请导出该状态方程。 题图 2-1 2i 和3u 存在标准状态方程 323 3212222))2(tan (231 dt u i dt du u i u i i di s =--+=- 2-2 题图2-2所示电路，非线性电阻的特性为：2 2223R R R u u i -=，试导出电路的状态方程。 题图 2-2 L C C L C C L C L s C i L R u L u L dt di u u C i C du i C i C du 2212 222221 111 1)3(11dt 1 1dt --=--=-= 2-3 试确定下列函数是否满足全局Lipschitz 条件 （1）2 211212()[2]T f x x x x x x =--可能不满足 （2）2 2 2 112()[]x x T f x x e x e --=满足

2-4 Van der pol 方程可以用状态方程描述为 122 2112(1)x x x x x x ε=??=-+-? 试证明，任取初始条件1020x x ，，对于某些充分小的δ，状态方程在[0]δ上有唯一解。 2-5 考虑标量微分方程 0tan(()),(0)x x t x x == 试证明微分方程对于任意0x ，在区间[0,)∞上具有唯一解。 2-6 已知非线性系统的状态方程为 ? ?????????-----=???? ??????-t te x x x x x t dt dx dt dx 22212131 213tanh 43 试判断该状态方程是否有唯一解。 当00,0t t t ≥>时有唯一解 2-7 试求下列电路状态方程的平衡点。 （1）???????+-=-=dxy by dt dy cxy ax dt dx （0，0） （2）???????+-=++-=222 2y x y x dt dy y x x y dt dx （0，0） （3）???????-==3x x dt dy y dt dx （0，0）；（1，0）；（-1，0） （4）???????-==1sin 2x dt dy y dt dx ,2,1,0) ,1();k 1±±=-k k ππ，（ （5）???????+-=-=+3 1dy by dt dy e dt dx y x (0,0)；0,0d b )d b ,d b (); d b ,d b (≠>- -d

第三章平面问题的有限元法作业及答案

第三章 平面问题的有限元法作业 1. 图示一个等腰三角形单元及其节点编码情况，设μ=0，单元厚度为t 。求 1）形函数矩阵[]N ；2）应变矩阵[]B ；3）应力矩阵[]S 。 4 第1题图 第2题图 2. 如题图所示，结构为边长等于a 的正方形，已知其节点位移分别为：11(,)u v 、 22(,)u v 、33(,)u v 、44(,)u v 。试求A 、B 、C 三点的位移。其中A 为正方形形心，B 为三角形形心。 3.直角边边长为l 的三角形单元，如题图所示。试计算单元等效节点载荷列阵（单元厚度为t ，不计自重）。 第3题图 第4题图 4. 如题图所示，各单元均为直角边边长等于l 的直角三角形。试计算（1）单元等效节点载荷列阵；（2）整体等效节点载荷列阵。已知单元厚度为t ，不计自重。

5．下列3个有限元模型网格，哪种节点编号更合理？为什么？ 9 34 6 7912 11 34 6 12142 (a) (b) (c) 第5题图 6.将图示结构画出有限元模型；标出单元号和节点号；给出位移边界条件；并计算半带宽（结构厚度为t ）。 2a (a) (b) 无限长圆筒 (c) 第6题图 7. 结构如图所示，已知结构材料常数E 和 ，单元厚度为t 。利用结构的对称性，采用一个单元，分别计算节点位移和单元应力。 第7题图

答案： 1. 1）形函数 i x N a = ， j y N a = ， 1m x y N a a =-- 2）应变矩阵 []1000101 000101011011B a -????=-??--???? 3）应力矩阵 []100010100 01 0111 110022 2 2S a ? ???-? ?=-????- -? ?? ? 2. A 点的位移为 ()2312A u u u = + ， ()231 2A v v v =+ B 点的位移为 ()24313B u u u u = ++ ， ()2431 3B v v v v =++ C 点的位移为 ()1223C a u u u = + ， ()C 1223 a v v v =+ 3. 单元等效节点载荷列阵为 {}11 11 00003 663 T e i j i j R q q q q ?? =++?? ?? 4. （2）整体等效节点载荷向量为 {}111100006 322T R qlt P qlt P P qlt qlt ?? =-???? 7. （1） 减缩后的整体刚度方程 22 12 2 1222 22221110222021102(1)2 2102x x b b ab R b ab b P v Et ab a b ab ab R v b a μμμ μμμμμμ---??- - ??????????--?????? -??? ?=????---+ +? ???? ?????????-????+?? ? ? 节点位移