曲悬臂梁网格细化

本示例说明如何使用*REFINE MESH关键字卡根据某个解变量细化四面体网格。结构是一个曲悬臂梁（图53），使用C3D10单元非常粗糙地网格化。梁的左侧在z方向完全固定，左下节点还在x和y方向固定，右下节点在y方向固定。9个力单位的载荷施加到梁右面节点上的+y方向。这导致z方向的法向应力如图所示。梁经历弯曲，在底部产生拉应力，在顶部产生压应力。梁的输入卡（circ10p.inp）是CalculiX测试套件的一部分。

**图53:** 粗网格中z方向的法向应力
$\begin{figure}\epsfig{file=opt1.eps,width=10cm}\end{figure}$

**图54:** 粗网格的误差估计器
$\begin{figure}\epsfig{file=opt2.eps,width=10cm}\end{figure}$

这里只复制输入卡中的step信息：

*STEP
*STATIC
*CLOAD
LOAD,2,1.
*NODE PRINT,NSET=FIX,TOTALS=ONLY
RF
*SECTION PRINT,SURFACE=Sfix,NAME=SP1
SOF,SOM
*NODE FILE
U
*EL FILE
S
*REFINE MESH,LIMIT=50.
S
*END STEP

它说明了几种获取反力的方法。一种方法是使用*NODE PRINT关键字卡请求在.dat文件中存储RF。它作用于节点集，在这种情况下是梁左表面上的所有节点。参数TOTALS=ONLY表示只打印力的总和，而不是各个分量。如果节点集相邻的单元没有分布载荷（ surface distributed loads (pressure) nor volumetric distribute loads (gravity, centrifugal forces)），则*NODE PRINT选项效果很好。否则，打印的RF值将包含这些力的一部分。

第二种方法是定义一个面部表面，并在*SECTION PRINT卡下使用SOF和SOM来请求该表面的力和力矩。表面Sfix由梁左表面上所有面组成。力和力矩通过跨表面积分获得。

.dat文件中的输出如下：

 total force (fx,fy,fz) for set FIX and time  0.1000000E+01

       -9.372063E-13 -9.000000E+00  3.127276E-12


 statistics for surface set SFIX and time  0.1000000E+01

   total surface force (fx,fy,fz) and moment about the origin(mx,my,mz)

    2.454956E+00 -7.226251E+00  1.377949E+01  7.236961E+01 -5.740438E+00 -4.957194E+00

   center of gravity and mean normal

    5.000000E-01  5.000000E-01  0.000000E+00  0.000000E+00  0.000000E+00 -1.000000E+00

   moment about the center of gravity(mx,my,mz)

    6.547987E+01  1.149306E+00 -1.165902E-01

   area,  normal force (+ = tension) and shear force (size)

    6.000000E+00 -1.377949E+01  7.631875E+00

由此可以观察到，*NODE PRINT语句获得的反力非常准确，然而跨表面应力的积分相当不准确：得到7.23个力单位而不是9个。绕重心的力矩是65.5 [力][长度]，而预期是72 [力][长度]（梁的长度为8个长度单位）。

误差估计器的值如图54所示。不出所料，误差相当高，高达30%。

为了获得更好的结果，通过*REFINE MESH,LIMIT=50卡触发自动基于应力的细化。细化的字段列在此卡下面。"S"表示将使用米塞斯应力。对于本示例，米塞斯应力达到约400应力单位，因此局部细化因子最高可达8（选择了50的细化限制）。在当前版本的CalculiX中，最多执行三次迭代，每次允许细化因子为2。细化总是应用于原始网格的版本，其中任何二次单元都被线性单元替换（C3D10由C3D4），即不考虑中间节点。这些细化迭代的结果存储为输入卡（仅包含网格）在文件finemesh.inp0、finemesh.inp1和finemesh.inp2中。生成存储在finemesh.inp2中的网格后，如果输入卡至少包含一个二次单元，程序将为所有单元生成中节点。所有节点随后投影到原始网格的面上。这意味着几何形状基本上由输入卡中网格的外表面描述。除四面体单元外的其他单元保持不变。结果的投影网格作为输入卡存储在jobname.fin中。它只包含细化的网格（节点和单元）。

运行circ10p输入卡并重新应用必要的边界和加载条件（必须手动完成）导致输入卡cric10pfin.inp（也是CalculiX测试示例的一部分）。运行此卡得到图55中的z方向法向应力和图56中的误差。

**图55:** 细网格中z方向的法向应力
$\begin{figure}\epsfig{file=refine3.eps,width=10cm}\end{figure}$

**图56:** 细网格的误差估计器
$\begin{figure}\epsfig{file=refine4.eps,width=10cm}\end{figure}$

网格在梁左侧面附近（应力最高的地方）进行了细化。结果的单元是二次单元，原始网格的曲率很好地保持了。

压应力略有增加，而拉应力现在更多地集中在y方向固定的节点附近。然而，总体水平相似。应力误差与粗网格大致相同，然而，在应力高的位置，误差现在较低，约为5%而不是30%。这些是感兴趣的位置。

.dat文件中反力的输出如下：

 total force (fx,fy,fz) for set FIX and time  0.1000000E+01

        3.221013E-12 -9.000000E+00  7.356782E-12


 statistics for surface set SFIX and time  0.1000000E+01

   total surface force (fx,fy,fz) and moment about the origin(mx,my,mz)

    1.512388E-01 -9.252627E+00 -7.227514E-01  7.175724E+01  1.563390E-01 -4.206416E+00

   center of gravity and mean normal

    5.000000E-01  5.000000E-01  4.014218E-19 -4.263022E-20  4.286885E-20 -1.000000E+00

   moment about the center of gravity(mx,my,mz)

    7.211862E+01 -2.050367E-01  4.955169E-01

   area,  normal force (+ = tension) and shear force (size)

    6.000000E+00  7.227514E-01  9.253863E+00

节点输出 again非常准确，而截面输出明显改善：总反力现在是-9.25个力单位，绕重心的力矩是72.12 [力][长度]。更细的网格导致更准确的节点应力，而这些正是用于确定截面力的。