数字旗手

电气化、自动化、数字化、智能化、智慧化

0%

引子 这个例子的亮点:介绍了最简单的Laplace方程的组装矩阵和右端项的一键生成函数。 这个例子来求解仅考虑Neumann边界条件的Laplace问题: 如果有解的话,解必须满足协调条件: 这里考虑计算域是以原点为圆心,半径为1的圆,$f=-2,g=1$是满足协调条件的。 虽然协调条件允许有解,但解仍然是不定的:在方程中仅仅解的导数固定,解可以是任意常数。所以需要施加另外一个条件来固定这个常数。可以有以下可能的方法: * 固定离散后的某个节点值为0或其他固定值。这意味着一个额外的条件$u_h(x_0)=0$。尽管这是通常的做法,但不是一个好方法,因为我们知道Laplace方程的解是在
Read more »

引子 本例将要完成以下目标: * 求解对流方程$\beta \cdot \nabla u = f$ * 使用多线程求解 * 设计一个简单的细化准则 方程离散 对流方程中$\beta$是一个描述对流方向和速度的矢量场,可能与空间位置相关,$f$是源项,$u$是解。这个方程求解的物理过程可以是浓度u在给定流场下的传输,没有扩散过程,但有源。 明显地,在入流侧,上述方程需要在边界上获得补充: 其中,$\Omega_-$代表边界的入流部分,其定义为: 其中,${\mathbf n}({\mathbf x})$是点$x$上的外法线。这种定义方式是非常直观的:因为如果外法线是朝外的,那么在入
Read more »

引子 真实生活中,大多数的偏微分方程都是一组方程,相应地,解也通常是矢量场。跟之前单方程的求解以及解是标量场相比,这种问题只是在组装矩阵和右端项时有些复杂,其他都一样。 本例中求解的是弹性问题: 其中,$c_{ijkl}$是刚度系数,通常与空间坐标相关。 弹性方程是Laplace方程的一种扩散形式,其解$u_l$是矢量场,表示一个刚体在力的作用下的位移。力$f_i$也是矢量场。 $c_{ijkl}$是一个四阶张量,共有$3^4=81$个分量,但其实只有21个分量是独立的,这是因为: 首先因为$\sigma_{ij}$和$\epsilon_{ij}$都是对称张量,导致$c_{ijkl}=c_{
Read more »

引子 在本例中,将会着眼于以下两方面: 1. 验证程序的正确性,生成收敛性统计表格; 2. 对于Helmholtz方程施加非齐次Neumann边界条件。 另外还有一些小的优化点。 验证程序正确性 也许从来不会有任何一个有限元程序一开始就是正确的,所以找到方法来验证计算的解是否正确就很有必要。通常选择已知精确解析解,并且比较精确解析解和计算离散解两者之间差别来求证。如果随着误差次数提高,两者之间差别逐渐趋于0,就说明程序的正确性。deal.II中就提供了这样一个函数:VectorTools::integrate_difference(),它提供了很多种范数的计算: 这些公式也适用于矢
Read more »

引子 本例主要着眼于处理局部细化的网格。如果临近单元细化多次以后,单元界面上的格点可能在另一边就不平衡,称为“悬点”。为了保证全局解在这些格点上也是连续的,必须对这些节点上的解的值施加一些限制,相应的核心的类是ConstraintMatrix。 程序解析 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 #include #include #include
Read more »

引子 此例没有介绍革命性的功能,但有很多对前面例子的“微创新”,包括: * 在不断细化的网格上的计算。数值计算通常要在不同的网格上进行,这样才能感受到精度。而且deal.II支持自适应网格,虽然这个例子中没有用到,但基础在这 * 读入非规则网格数据 * 计算优化 * debug模式,使用Assert宏 * 变系数Possion方程,使用预条件迭代求解器 这里要求解的方程是: 如果$a(x)$是常系数,那么就成了Possion方程,如果它是空间相关的变系数,方程就复杂一些了。 还是得先写出方程的弱形式: 程序解析 以下是头文件们: 1 2 3 4 5 6 7 8 9 10 11
Read more »

引子 deal.II有一个特性,叫作”维度无关的编程“。前面的例子中,很多类都有一个尖括号括起的数字的后缀。 这意味着该类能作用在不同的维度上,而不同维度的计算代码基本相同,这能显著地减少重复代码。这正是C++的模板template的拿手好戏。 在Step4中,将显示程序怎样维度无关的编程,具体是将step3中的Laplace问题同时在二维和三维上求解,以及右端项不再是常量、边界值不再为0。 程序解析 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 #include
Read more »

引子 这是使用有限元法进行具体计算的第一个算例,求解的是一个简化的Possion方程,其在边界上为0,而右端项不为0,即: 求解域是单位正方形$\Omega=[0,1]^2$,其上的网格划分在step1和step2中已经涉及。 这里也仅仅计算特例$f(x)=1$,更一般的情形详见step4。 推导 首先需要得到上述方程的弱形式,即在方程两侧左乘一个测试函数$\phi$并在计算域上积分(注意是左乘,而不是右乘): 然后利用高斯散度公式进行分部积分,可得: 同时,因为测试函数必须满足相同的边界条件,即在边界上$\phi=0$,所以上式变为: 这里应用了通常的简写形式:$(a,b)=\in
Read more »

本文是对Mathematica有限元分析与工程应用一书的学习笔记。 三角形单元适用于具有复杂和不规则边界形状的问题,是一种常见的网格离散方式。 双线性三角形单元 局部坐标系 之前的杆单元和桁架元的局部坐标系容易建立,即建在其自身上即可。而三角形单元的局部坐标系显然不能这样建立,其常采用一套无量纲的自然坐标系——面积坐标,如下图所示: 三角形123内部有一任意点P,P与顶点1、2、3组成3个子三角形,每个子三角形的面积与总面积之比记为$L_i$,即P点的面积坐标为$(L_1,L_2,L_3)$。 因为三个坐标相加为1,所以仅有两个独立的自然坐标,所以可以简记为: 注意到:面积坐标还有如下
Read more »

本文是对Mathematica有限元分析与工程应用一书的学习笔记。 桁架元的特点 平面桁架元是既有局部坐标又有整体坐标的二维有限元,因此比起之前的杆单元,需要多一步坐标变换。 桁架元示意图如下: 指定整体坐标系为X-Y,局部坐标系为x-y。则两者之间的转换关系为: 即: 其中: 局部坐标系下的有限元方程为: 为了把有限元方程从局部坐标系变换到整体坐标系,可通过转换矩阵: 所以: 又因为转换矩阵T满足如下关系(可实际计算验证一下): 所以: 所以整体坐标系的刚度矩阵与局部坐标系的刚度矩阵关系为: 由于局部坐标系下的单元刚度矩阵为: 其中$k=\frac{EA}{L}$。那
Read more »