求解偏微分方程开源有限元软件deal.II学习--Step 8

最新推荐文章于 2024-07-12 15:51:28 发布

北冥有鱼虚空充盈

最新推荐文章于 2024-07-12 15:51:28 发布

阅读量1k

点赞数

分类专栏： finite element method

finite element method 专栏收录该内容

34 篇文章 17 订阅

订阅专栏

求解偏微分方程开源有限元软件deal.II学习--Step 8

Posted on 2016-09-06 | In computational material science | 4条评论

引子

真实生活中，大多数的偏微分方程都是一组方程，相应地，解也通常是矢量场。跟之前单方程的求解以及解是标量场相比，这种问题只是在组装矩阵和右端项时有些复杂，其他都一样。
本例中求解的是弹性问题：

- \partial j (c i j k l \partial k u l) = f i, i = 1 \dots d

其中， cijkl 是刚度系数，通常与空间坐标相关。
弹性方程是Laplace方程的一种扩散形式，其解 ul 是矢量场，表示一个刚体在力的作用下的位移。力 fi 也是矢量场。
cijkl 是一个四阶张量，共有 34=81 个分量，但其实只有21个分量是独立的，这是因为：
首先因为 σij 和 ϵij 都是对称张量，导致 cijkl=cjikl 和 cijkl=cijlk ，这样就将独立的弹性常数减少到 (3∗2∗1)2=36 个。
然后因为应变能密度函数与应力应变的关系，导致 cijkl=cklij ，这样代表沿对角线对称，这样将弹性常数减少到21个。
这样三维情形下应力应变关系矩阵形式为：

⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ σ 11 σ 22 σ 33 σ 12 σ 23 σ 31 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ c 1111 c 1122 c 2222 c 1133 c 2233 c 3333 c 1112 c 2212 c 3312 c 1212 c 1123 c 2223 c 3323 c 1223 c 2323 c 1131 c 2231 c 3331 c 1231 c 2331 c 3331 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ϵ 11 ϵ 22 ϵ 33 2 ϵ 12 2 ϵ 23 2 ϵ 31 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥

在各向同性材料中，通过引入两个系数，系数张量变成(可以通过理论证明，各向同性的均匀弹性体的弹性常数只有两个)：

c i j k l = λ δ i j δ k l + μ (δ i k δ j l + δ i l δ j k) .

比如，当 i=j=k=l=1 时：

c 1111 = λ + 2 μ

当 i=j=k=1,l=2 时，

c 1112 = 0

当 i=k=1,j=l=2 时，

c 1212 = μ

当 i=j=1,k=l=2 时，

c 1122 = λ

这就是各向同性材料的应力应变关系矩阵：

⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ λ + 2 μ λ λ + 2 μ λ λ λ + 2 μ 000 μ 0000 μ 00000 μ ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥

回到力与位移的系数矩阵，按各个字母的序号循环可得：

它的系数矩阵是一个 2∗2 的矩阵，其中各个元素可由上图得到。
对于各向同性材料，代入上面的取值后，有：

- \nabla λ (\nabla \cdot u) - (\nabla \cdot μ \nabla) u - \nabla \cdot μ (\nabla u) T = f,

其展开就是如图：

注意各种标量、矢量和张量的梯度和散度运算。
其弱形式为：

下面就是怎样组装这个线性系统。第一件事情就是需要知道在矢量值的有限元中形函数是怎样工作的。大体过程这样：设 n

为要建立的矢量单元的分量，即标量单元，的形函数的个数，比如之前用的双线性单元，二维情形下 n=4 。设 N 是矢量单元的形函数的个数，二维情形下， N=2n ，那么矢量单元的第i个形函数的形式为：

Φ i (x) = φ base (i) (x) e comp (i),

其中， comp(i) 是告诉我们 Φi 的哪个分量非0(对每一个矢量形函数，只有一个分量非0，其余分量都为0)，比如 comp(1)=0 表示第1个矢量形函数的第0个分量非0。 φbase(i)(x) 描述形函数与坐标的关系，就是标量单元的第 base(i) 个形函数，比如具体形式是这样的：

Φ 0 (x) Φ 1 (x) Φ 2 (x) Φ 3 (x) = [φ 0 (x) 0], = [0 φ 0 (x)], = [φ 1 (x) 0], = [0 φ 1 (x)], \dots

其中：

comp (0) = 0, comp (1) = 1, comp (2) = 0, comp (3) = 1, \dots base (0) = 0, base (1) = 0, base (2) = 1, base (3) = 1, \dots

在绝大多数情况下，不需要知道哪个 φbase(i) 属于 Φi ，于是定义：

ϕ i = φ base (i)

所以，矢量形函数表示为：

Φ i (x) = ϕ i (x) e comp (i) .

使用上述矢量形函数，构造出离散有限元解：

u h (x) = \sum i Φ i (x) U i

其中， Ui 是系数，是标量。定义一个类似的函数 vh 作为试探函数，那么问题就变为：找到系数 Ui ，使得：

a (u h, v h) = (f, v h) \forall v h .

将双线性的具体形式代入，可得：

\sum i, j U i V j \sum k, l [(λ \partial l (Φ i) l, \partial k (Φ j) k) Ω + (μ \partial l (Φ i) k, \partial l (Φ j) k) Ω + (μ \partial l (Φ i) k, \partial k (Φ j) l) Ω] = \sum j V j \sum l (f l, (Φ j) l) Ω .

注意到：下标k和l是对所有空间方向进行循环， 0≤k,l<d ，而下标i和j是对所有自由度进行循环。
那么，单元K上的单元刚度矩阵就是：

A K i j = \sum k, l [(λ \partial l (Φ i) l, \partial k (Φ j) k) K + (μ \partial l (Φ i) k, \partial l (Φ j) k) K + (μ \partial l (Φ i) k, \partial k (Φ j) l) K]

这里，i和j是局部自由度，有 0≤i,j<N 。
在这些公式中，我们通常取矢量形函数的部分分量，根据定义，有：

(Φ i) l = ϕ i δ l, comp (i),

那么，进一步简化得到：

A K i j = \sum k, l [(λ \partial l ϕ i δ l, comp (i), \partial k ϕ j δ k, comp (j)) K + (μ \partial l ϕ i δ k, comp (i), \partial l ϕ j δ k, comp (j)) K + (μ \partial l ϕ i δ k, comp (i), \partial k ϕ j δ l, comp (j)) K] = (λ \partial comp (i) ϕ i, \partial comp (j) ϕ j) K + \sum l (μ \partial l ϕ i, \partial l ϕ j) K δ comp (i), comp (j) + (μ \partial comp (j) ϕ i, \partial comp (i) ϕ j) K = (λ \partial comp (i) ϕ i, \partial comp (j) ϕ j) K + (μ \nabla ϕ i, \nabla ϕ j) K δ comp (i), comp (j) + (μ \partial comp (j) ϕ i, \partial comp (i) ϕ j) K .

同样地，单元对右端项的贡献为：

f K j = \sum l (f l, (Φ j) l) K = \sum l (f l, ϕ j δ l, comp (j)) K = (f comp (j), ϕ j) K .

程序解析

#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>

以上是以前用过的头文件。

1	#include <deal.II/fe/fe_system.h>

该头文件提供对矢量值的有限元的支持。

1	#include <deal.II/fe/fe_q.h>

从常规的Q1单元组合得到矢量值的有限元，Q1单元在以上头文件中。

1 2	#include <fstream> #include <iostream>

C++的标准库。
建立step8的命名空间：

namespace Step8
{
    using namespace dealii;
template <int dim>
    class ElasticProblem
    {
        public:
            ElasticProblem ();
            ~ElasticProblem ();
            void run ();
        private:
            void setup_system ();
            void assemble_system ();
            void solve ();
            void refine_grid ();
            void output_results (const unsigned int cycle) const;
            Triangulation<dim> triangulation;
            DoFHandler<dim> dof_handler;
            FESystem<dim> fe;
            ConstraintMatrix hanging_node_constraints;
            SparsityPattern sparsity_pattern;
            SparseMatrix<double> system_matrix;
            Vector<double> solution;
            Vector<double> system_rhs;
    };

step8的类跟step6差不多，唯一一个变化是fe的类型，这里使用的是FESystem，不再是FE_Q。实际上FESystem本身不是一个有限元类型，不提供形函数。它就是将多个单元集合起来形成一个矢量的有限单元。
然后建立右端项：

template <int dim>
class RightHandSide : public Function<dim>
{
    public:
        RightHandSide ();
        virtual void vector_value (const Point<dim> &p,
                Vector<double> &values) const;
        virtual void vector_value_list (const std::vector<Point<dim> > &points,
                std::vector<Vector<double> > &value_list) const;
};

vector_value是取得某个位置的矢量值，vector_value_list是一下取得很多。

template <int dim>
RightHandSide<dim>::RightHandSide ()
    :
        Function<dim> (dim)
{}

析构函数中给Function传递了参数，代表分量的个数，这里就是dim。

template <int dim>
inline
void RightHandSide<dim>::vector_value (const Point<dim> &p,
        Vector<double> &values) const
{
    Assert (values.size() == dim,
            ExcDimensionMismatch (values.size(), dim));
    Assert (dim >= 2, ExcNotImplemented());

这里加入了几个提前判断，用来保证维数和矢量大小正确。

Point<dim> point_1, point_2;
point_1(0) = 0.5;
point_2(0) = -0.5;
if (((p-point_1).norm_square() < 0.2*0.2) ||
        ((p-point_2).norm_square() < 0.2*0.2))
    values(0) = 1;
else
    values(0) = 0;

如果在离两个圆心一定范围内，就把x方向的力设为1，否则设为0。

if (p.norm_square() < 0.2*0.2)
    values(1) = 1;
else
    values(1) = 0;
}

如果离原点一定范围内，就把y方向的力设为1，否则设为0。

template <int dim>
void RightHandSide<dim>::vector_value_list (const std::vector<Point<dim> > &points,
        std::vector<Vector<double> > &value_list) const
{
    Assert (value_list.size() == points.size(),
            ExcDimensionMismatch (value_list.size(), points.size()));
    const unsigned int n_points = points.size();
    for (unsigned int p=0; p<n_points; ++p)
        RightHandSide<dim>::vector_value (points[p],
                value_list[p]);
}

然后就是一下取得好多点上的值。

template <int dim>
ElasticProblem<dim>::ElasticProblem ()
    :
        dof_handler (triangulation),
        fe (FE_Q<dim>(1), dim)
{}
template <int dim>
ElasticProblem<dim>::~ElasticProblem ()
{
    dof_handler.clear ();
}

这是求解类的构造和析构函数。在构造函数中，给fe传递两个参数：一个是构造矢量有限元所基于的标量有限元，另一个就是多少个标量堆起来等于1个矢量，这里就是dim。知道这些信息后，FESystem就知道该怎么合成矢量有限元了。

template <int dim>
void ElasticProblem<dim>::setup_system ()
{
    dof_handler.distribute_dofs (fe);
    hanging_node_constraints.clear ();
    DoFTools::make_hanging_node_constraints (dof_handler,
            hanging_node_constraints);
    hanging_node_constraints.close ();
    DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
    DoFTools::make_sparsity_pattern(dof_handler,
            dsp,
            hanging_node_constraints,
            / *keep_constrained_dofs = * / true);
    sparsity_pattern.copy_from (dsp);
    system_matrix.reinit (sparsity_pattern);
    solution.reinit (dof_handler.n_dofs());
    system_rhs.reinit (dof_handler.n_dofs());
}

建立系统。过程跟step6相同，这里的这些类都能处理矢量的单元，实际无论矢量还是标量单元，这些类都一视同仁。
重头戏就是组装系统了。

template <int dim>
void ElasticProblem<dim>::assemble_system ()
{
    QGauss<dim> quadrature_formula(2);
    FEValues<dim> fe_values (fe, quadrature_formula,
            update_values | update_gradients |
            update_quadrature_points | update_JxW_values);
    const unsigned int dofs_per_cell = fe.dofs_per_cell;
    const unsigned int n_q_points = quadrature_formula.size();
    FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
    Vector<double> cell_rhs (dofs_per_cell);
    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);

第一部分跟之前的相同：设置合适的积分公式、初始化FEValues对象、声明一些附加数组等。其中，获取每个单元上的自由度数，是从合成过的有限单元上取得，而不是那个标量Q1对象。这里，自由度数等于dim乘以Q1单元的每个单元上的自由度数。

1 2	std::vector<double> lambda_values (n_q_points); std::vector<double> mu_values (n_q_points);

然后存储所有积分点上两个系数的值。

1	ConstantFunction<dim> lambda(1.), mu(1.);

这里将两个系数都设为定值。

1
2
3

RightHandSide<dim> right_hand_side;
std::vector<Vector<double> > rhs_values (n_q_points,
        Vector<double>(dim));

建立右端项的对象。因为是矢量，所以rhs_values的类型也变化了。
然后开始单元循环：

typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
         endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
    cell_matrix = 0;
    cell_rhs = 0;
    fe_values.reinit (cell);
    lambda.value_list (fe_values.get_quadrature_points(), lambda_values);
    mu.value_list (fe_values.get_quadrature_points(), mu_values);
    right_hand_side.vector_value_list (fe_values.get_quadrature_points(),
            rhs_values);

上面计算系数和右端项在积分点上的值。
然后就是计算单元刚度矩阵和右端项：

for (unsigned int i=0; i<dofs_per_cell; ++i)
{
    const unsigned int
        component_i = fe.system_to_component_index(i).first;
    for (unsigned int j=0; j<dofs_per_cell; ++j)
    {
        const unsigned int
            component_j = fe.system_to_component_index(j).first;
        for (unsigned int q_point=0; q_point<n_q_points;
                ++q_point)
        {
            cell_matrix(i,j)
                +=

整个计算过程完全对应前面引子中的推导。component_i就是非0分量的指标，它通过fe.system_to_component_index(i).first这个函数取得，实际上first取得矢量形函数的非0分量的指标，而second取得具体这个形函数的值，即引子中的base。

(
 (fe_values.shape_grad(i,q_point)[component_i] *
  fe_values.shape_grad(j,q_point)[component_j] *
  lambda_values[q_point])
 +
 (fe_values.shape_grad(i,q_point)[component_j] *
  fe_values.shape_grad(j,q_point)[component_i] *
  mu_values[q_point])
 +

这一项计算的是(lambda d_i u_i, d_j v_j) + (mu d_i u_j, d_j v_i)。

((component_i == component_j) ?
 (fe_values.shape_grad(i,q_point) *
  fe_values.shape_grad(j,q_point) *
  mu_values[q_point]) :
 0)
)
*
fe_values.JxW(q_point);
        }
    }
}

这一项计算的是(mu nabla u_i, nabla v_j)。注意这里的grad没有加后面的方括号，用了重载的乘号。
且使用了条件表达式，判断两个下标是否相同。
单元右端项的计算：

for (unsigned int i=0; i<dofs_per_cell; ++i)
{
    const unsigned int
        component_i = fe.system_to_component_index(i).first;
    for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
        cell_rhs(i) += fe_values.shape_value(i,q_point) *
            rhs_values[q_point](component_i) *
            fe_values.JxW(q_point);
}

组装：

cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
    for (unsigned int j=0; j<dofs_per_cell; ++j)
        system_matrix.add (local_dof_indices[i],
                local_dof_indices[j],
                cell_matrix(i,j));
    system_rhs(local_dof_indices[i]) += cell_rhs(i);
}
}
hanging_node_constraints.condense (system_matrix);
hanging_node_constraints.condense (system_rhs);

施加边界条件：

std::map<types::global_dof_index,double> boundary_values;
VectorTools::interpolate_boundary_values (dof_handler,
        0,
        ZeroFunction<dim>(dim),
        boundary_values);
MatrixTools::apply_boundary_values (boundary_values,
        system_matrix,
        solution,
        system_rhs);
}

这里做了一些小修改，因为解是矢量的，所以边界条件施加的也应该是矢量的。而ZeroFunction可以接收参数来形成不同类型的量，这里传递的是dim。
求解器：

template <int dim>
void ElasticProblem<dim>::solve ()
{
    SolverControl solver_control (1000, 1e-12);
    SolverCG<> cg (solver_control);
    PreconditionSSOR<> preconditioner;
    preconditioner.initialize(system_matrix, 1.2);
    cg.solve (system_matrix, solution, system_rhs,
            preconditioner);
    hanging_node_constraints.distribute (solution);
}

求解过程不变，求解器不管具体问题是什么，只要线性系统是正定且对称的，CG算法就能用。
细化网格：

template <int dim>
void ElasticProblem<dim>::refine_grid ()
{
    Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
    KellyErrorEstimator<dim>::estimate (dof_handler,
            QGauss<dim-1>(2),
            typename FunctionMap<dim>::type(),
            solution,
            estimated_error_per_cell);
    GridRefinement::refine_and_coarsen_fixed_number (triangulation,
            estimated_error_per_cell,
            0.3, 0.03);
    triangulation.execute_coarsening_and_refinement ();
}

指示子是用的所有方向上的位移具有相同的权重。
结果输出：

template <int dim>
void ElasticProblem<dim>::output_results (const unsigned int cycle) const
{
    std::string filename = "solution-";
    filename += ('0' + cycle);
    Assert (cycle < 10, ExcInternalError());
    filename += ".vtk";
    std::ofstream output (filename.c_str());
    DataOut<dim> data_out;
    data_out.attach_dof_handler (dof_handler);
    std::vector<std::string> solution_names;
    switch (dim)
    {
        case 1:
            solution_names.push_back ("displacement");
            break;
        case 2:
            solution_names.push_back ("x_displacement");
            solution_names.push_back ("y_displacement");
            break;
        case 3:
            solution_names.push_back ("x_displacement");
            solution_names.push_back ("y_displacement");
            solution_names.push_back ("z_displacement");
            break;
        default:
            Assert (false, ExcNotImplemented());
    }
    data_out.add_data_vector (solution, solution_names);
    data_out.build_patches ();
    data_out.write_vtk (output);
}

因为结果是矢量，所以给每个分量都有一个名字。
运行函数：

template <int dim>
void ElasticProblem<dim>::run ()
{
    for (unsigned int cycle=0; cycle<8; ++cycle)
    {
        std::cout << "Cycle " << cycle << ':' << std::endl;
        if (cycle == 0)
        {
            GridGenerator::hyper_cube (triangulation, -1, 1);
            triangulation.refine_global (2);
        }
        else
            refine_grid ();
        std::cout << " Number of active cells: "
            << triangulation.n_active_cells()
            << std::endl;
        setup_system ();
        std::cout << " Number of degrees of freedom: "
            << dof_handler.n_dofs()
            << std::endl;
        assemble_system ();
        solve ();
        output_results (cycle);
    }
}
}

Attention!!
这里有个小问题，刚开始产生网格后，就全局细化了两次。这是因为这里选择的右端项相当局域化，如果只细化一次，网格的积分点很稀疏，它捕捉的右端项的值全是0，这样就计算错误了。所以，要考虑到网格细化对值的正确捕捉。
Attention完毕！！
main函数如下：

int main ()
{
    try
    {
        Step8::ElasticProblem<2> elastic_problem_2d;
        elastic_problem_2d.run ();
    }
    catch (std::exception &exc)
    {
        std::cerr << std::endl << std::endl
            << "----------------------------------------------------"
            << std::endl;
        std::cerr << "Exception on processing: " << std::endl
            << exc.what() << std::endl
            << "Aborting!" << std::endl
            << "----------------------------------------------------"
            << std::endl;
        return 1;
    }
    catch (...)
    {
        std::cerr << std::endl << std::endl
            << "----------------------------------------------------"
            << std::endl;
        std::cerr << "Unknown exception!" << std::endl
            << "Aborting!" << std::endl
            << "----------------------------------------------------"
            << std::endl;
        return 1;
    }
    return 0;
}