水坝平面应变分析MATLAB与ANSYS对比（需要完整MATLAB代码的后台私信）
问题描述
等截面挡水坝如图1所示（梯形截面），坝底部长3m，顶部长1.2m，高5m。材料为混凝土，密度1930kg/m²，弹性模量E=30MPa，泊松比0.3，假设混凝土抗拉压作用相同坝体底部固定在地基上，斜面上受水压作用，水面距离坝底高4m。左下角顶点为坐标原点，计算水坝横截面应力分布，并给出x=2.5m垂直线上的x方向正应力分布。
图1 等截面挡水坝示意图
前处理
1.单位定义
ANSYS Mechanical APDL中无需指定单位，只需要保证物理量量纲统一即可。
2.定义单元类型
Main Menu中，Preprocessor -> Element Type -> Add/Edit/Delete，
在Element Types对话框中点击Add…，在Library of Element Types对话框中选择单元类型为Solid，Quad中的4 node 182，然后OK，关闭Library of Element Types对话框。
可见Element Types对话框中已经把该单元包含进来了，打开Options…选项，将K3选项从默认的Plane stress改为Plane strain，设为平面应变。OK关闭设置，Close关闭Element Types对话框。
3. 定义材料属性
Main Menu中，Preprocessor -> Material Props -> Material Models，
在打开的Define Material Model Behavior窗口中，选择Structural -> Linera -> Elastic -> Isotropic，在打开的对话框中的EX输入弹性模量为30e6，单位为Pa，PRXY为泊松比，指定为0.3，OK确认。
选择Structural -> Density，输入密度为1930，单位是kg/m3，OK确认。
4.创建模型
Main Menu中，Preprocessor -> Modeling -> Create -> Keypoints -> In Active CS，
在Create Keypoints in Active Coordinate System对话框中，
输入节点编号1，坐标0，0，0，Apply；
输入节点编号2，坐标3，0，0，Apply；
输入节点编号3，坐标3，5，0，Apply；
输入节点编号4，坐标1.8，5，0，OK。
Main Menu中，Preprocessor -> Modeling -> Create -> Areas -> Arbitirary -> Through KPs，
用鼠标依次选择1、2、3、4点，OK。
5.划分单元
Main Menu中，Preprocessor -> Meshing -> MeshTool，
在弹出的MeshTool对话框中，先确定网格尺度，点击Size Controls组类中Lines后面的Set，
在弹出的Element Size on Picked Lines对话框中，选择上下两条边，OK，
在弹出的Element Sizes on Picked Lines对话框中，设定单元划分数目NDIV为20，OK确定。
同样方法，将左右两条边的网格数目NDIV设置为40。
Main Menu中，Preprocessor -> Meshing -> MeshTool，
在MeshTool对话框中，保持Mesh：中下拉选项为Areas，保持Shape中的Quad，选择划分类型为Mapped，点击最下方的Mesh，在弹出的Mesh Lines点选框中，Pick All。
即可完成网格划分，如图2。
图2 网格划分示意图
三、求解
1.设置分析类型
Main Menu中，Solution -> Analysis Type -> New Analysis，保持默认的Static静态分析，OK。
2 施加约束（施加边界条件）
Main Menu中，Solution -> Define Loads -> Apply -> Structural -> Displacement -> On Lines，
在弹出的Apply U,ROT on Lines选择框中，用鼠标在图形窗口中点选最下部的水平直线段，OK，
在弹出的Apply U,ROT on Lines对话框中，DOFs to be constrained（限定自由度）选项中，选中All DOF，即限制该线上节点在X和Y方向的位移，OK确定。
3. 施加载荷
这里的载荷是随着水深变化的压力，所以要采用如下方法来添加它。
Utility Menu，Parameters -> Functions -> Define/Edit，弹出Function Editor窗口，
在数字键上方的下拉列表中选择Y作为自变量，则Results文本框中自动出现了{Y}；
在Results文本框中输入随着深度变化的压力函数9800*(4-{Y})，即p = ρ × g × h = 9800 × ( 4 − y )；
File->Save保存该函数，将其命名为“waterPress.func”，关闭Function Editor窗口。
Utility Menu，Parameters -> Functions -> Read From File，弹出的窗口中选择刚才保存的函数文件waterPress.func，打开。
在跳出的Function Loader窗口中，将其命名为“MyPress”即可，OK。
Main Menu中，Solution -> Define Loads -> Apply -> Structural -> Pressure -> On Nodes，
跳出的Apply Pres on Lines窗口中用鼠标拾取最左侧谢直线上从下到上的33个节点，OK。（注意，由于这条斜线均分了40个单元，共有41个节点，水面以下收到压力载荷作用，这条线总共高度是5m，水面深度是4m，所以水面以下的节点是从下到上共（40*4/5+1）=33个节点）
跳出的Apply Pres on Nodess窗口中，将Apply PRES on nodes as a选择为Existing table，OK。
Main Menu，Solution -> Define Loads -> Apply -> Structural -> Inertia -> Gravity -> Global，ACELY中输入9.8，OK。
Main Menu中，Solution -> Solve -> Current LS，跳出的/STATUS Command是对当前待求解问题的描述，直接点击Solve Current Load Step中的OK
四、后处理
1.查看变形状态
Main Menu中，General PostProc -> Plot Results -> Deformed Shape，
Main Menu中，General PostProc -> Plot Results -> Deformed Shape，
在Plot Deformed Shape窗口中，选择Def + undeformed，OK。
不难发现，在左下方水压的作用下，坝体截面的变形就恰如一端固支梁的弯曲问题，其整体向右侧发生了弯曲变形，且y方向略有压缩（重力作用会使得y方向有压缩，而x方向弯曲后y方向也会有压缩效应），如图3所示。
图3 水坝变形状态
2.绘制变形和应力云图
Main Menu中，General PostProc -> Plot Results -> Contour Plot -> Nodal Solu，
选择Stress中的X-Component of stress，OK
对于底部分析发现，底边前半部分受拉而后半部分受压，而最大应力位于右下角点附近，这符合物理直觉。即，一个弯曲的梁，在其固定端处的截面上，一定是有部分受压缩，而有部分受拉伸作用，如图4所示。
图4 水坝变形和应力云图
程序设计思路
本程序采用二维平面应变有限元方法，对梯形混凝土重力坝在自重与水压力作用下的应力分布进行求解。首先，根据坝体几何建立结构网格，并生成四节点矩形单元划分；然后基于材料参数构造平面应变条件下的本构关系矩阵，通过高斯积分计算每个单元的刚度矩阵并组装成整体刚度矩阵；接着施加水压力和重力载荷，并对底部及边界进行约束处理；最后求解线性方程组得到位移，再由位移计算各单元中心的x方向应力（σx），并绘制应力云图展示坝体内部的应力分布情况。
