一个简单的面向对象有限元程序的开发

发表于 | 分类于
\( \def\ <#1>{\left <#1\right>} \newcommand{\CC}{\bm{C}} \newcommand{\dydx}[2]{\frac{\mathrm{d}#1}{\mathrm{d}#2}} \newcommand{\pypx}[2]{\frac{\partial #1}{\partial #2}} \newcommand{\pyypxx}[2]{\frac{\partial^2 #1}{\partial #2^2}} \newcommand{\dyydxx}[2]{\frac{\mathrm{d}^2 #1}{\mathrm{d} #2^2}} \newcommand{\half}{\frac{1}{2}} \newcommand{\pprime}{\prime\prime} \) \( \newcommand{\bm}[1]{\boldsymbol{\mathbf{#1}}} \)

本程序使用C++编写，GUI界面使用Qt库，线性代数运算使用Eigen库，实现了二维桁架结构、钢架结构的通用有限元程序，能够计算静载荷下的位移，以及固有频率和固有振型。

理论基础

桁架结构

二维桁架结构的的刚度矩阵为

质量矩阵为

其中，$\rho$为线密度。坐标变换矩阵为

那么在总体坐标系下的质量矩阵，载荷向量，位移向量为 \begin{equation} \bm{\tilde m}=\bm{T}^T\bm{mT}=\bm{m},\quad \bm{\tilde f}=\bm{T}^T\bm{f},\quad \bm{\tilde v}=\bm{T}^T\bm{v} \end{equation}

弯曲单元

采用三次梁单元，建立梁单元的刚度矩阵为

单元质量矩阵为

对于二维平面上的梁，有坐标转换矩阵为

那么在总体坐标系下的刚度矩阵为 \begin{equation} \bm{\tilde k}=\bm{T}^T\bm{kT} \end{equation} 同理有 \begin{equation} \bm{\tilde m}=\bm{T}^T\bm{mT},\quad \bm{\tilde f}=\bm{T}^T\bm{f},\quad \bm{\tilde v}=\bm{T}^T\bm{v} \end{equation}

约束的处理

如果无约束的有限元方程是下面的形式：

那么应用约束条件后，有限元方程应该是这样的形式：

有限元部分代码

class trussElement;

class node
{
    friend class trussElement;
    friend class trussStructure;
public:
    node(double x1, double y1)
        :m_x(x1), m_y(y1) {}

private:
    double m_x, m_y;
};

struct Constraint
{
    bool UX;
    bool UY;
    int node;
    Constraint(bool X, bool Y, int n){
        UX = X; UY = Y;
        node = n;
    }
};

std::vector<node> globalNodeList;
// Here I should use singleton later

class trussElement
{
    EIGEN_MAKE_ALIGNED_OPERATOR_NEW
public:
    trussElement(std::vector<node>& nodeList,
                 int node1, int node2, double EA, double rho);

    int getNumNode1() {return mNumNode1;}
    int getNumNode2() {return mNumNode2;}

//private:
    int mNumNode1;
    int mNumNode2;

    double m_theta;
    double m_x1, m_x2, m_y1, m_y2;
    double m_EA, m_rho;
    double m_l;

    Eigen::Matrix<double, 4, 4> m_K;
    Eigen::Matrix<double, 4, 4> m_mass;
};

trussElement::trussElement(std::vector<node>& nodeList,
                           int node1, int node2, double EA, double rho) :
    m_x1(nodeList[node1].m_x), m_y1(nodeList[node1].m_y),
    m_x2(nodeList[node2].m_x), m_y2(nodeList[node2].m_y),
    mNumNode1(node1), mNumNode2(node2),
    m_EA(EA), m_rho(rho)
{
    m_theta = std::atan2(m_y2 - m_y1, m_x2 - m_x1);
    m_l = sqrt((m_y2 - m_y1) * (m_y2 - m_y1) 
            + (m_x2 - m_x1) * (m_x2 - m_x1));

    double c = cos(m_theta);
    double s = sin(m_theta);

    m_K << c * c, s * c, -c * c, -s * c,
           s * c, s * s, -s * c, -s * s,
           -c * c, -s * c, c * c, s * c,
           -s * c, -s * s, s * c, s * s;

    m_K *= m_EA / m_l;

    Eigen::Matrix<double, 4, 4> mass0, trans;
    mass0 << 2., .0, 1., .0,
             .0, 2., .0, 1.,
             1., .0, 2., .0,
             .0, 1., .0, 2.;
    mass0 *= m_rho * m_l / 6.0;

    trans << c, s, .0, .0,
             -s, c, .0, .0,
             .0, .0, c, s,
             .0, .0, -s, c;

    m_mass = trans.transpose() * mass0 * trans;
}

class trussStructure
{
public:
    EIGEN_MAKE_ALIGNED_OPERATOR_NEW
    void addTrussElement(const trussElement& elem);

    void addConstraint(Constraint con);

    void removeTrussElement(int i);

    void recalcAllMatrices();

    void addLoad(int numNode, double Fx, double Fy);

//private:
    void recalcLoad();
    void recalcStiffness();
    void recalcMass();

    std::vector<trussElement> trussArray;
    std::vector<Constraint> constrainsArray;
    std::map<int, double> xLoadOnNode;
    std::map<int, double> yLoadOnNode;
    std::vector<int> indicesToConstraint;

    // assemble matrices are no constraints applied,
    // while global matrices are applied.
    Eigen::MatrixXd m_assembleStiffness;
    Eigen::MatrixXd m_assembleMass;
    Eigen::MatrixXd m_globalStiffness;
    Eigen::MatrixXd m_globalMass;
    Eigen::MatrixXd m_assembleLoad;
    Eigen::MatrixXd m_globalLoad;


    Eigen::MatrixXd m_displacement;
};

void trussStructure::addTrussElement(const trussElement &elem)
{
    trussArray.push_back(elem);
}

void trussStructure::addConstraint(Constraint con)
{
    constrainsArray.push_back(con);
    if (con.UX){
        indicesToConstraint.push_back(2 * con.node + 0);
    }
    if (con.UY){
        indicesToConstraint.push_back(2 * con.node + 1);
    }
}

void trussStructure::addLoad(int numNode, double Fx, double Fy)
{
    xLoadOnNode[numNode] += Fx;
    yLoadOnNode[numNode] += Fy;
}

void trussStructure::recalcAllMatrices()
{
    std::sort(indicesToConstraint.begin(), indicesToConstraint.end());
    recalcLoad();
    recalcMass();
    recalcStiffness();
}

void trussStructure::recalcStiffness()
{
    int n = globalNodeList.size();
    m_assembleStiffness = Eigen::MatrixXd::Zero(2 * n, 2 * n);
    for(auto elem : trussArray){
        Eigen::MatrixXd mapMatrix = Eigen::MatrixXd::Zero(4, 2 * n);
        int a1 = elem.getNumNode1();
        int a2 = elem.getNumNode2();
        mapMatrix(0, 2 * a1) = 1.0;
        mapMatrix(1, 2 * a1 + 1) = 1.0;
        mapMatrix(2, 2 * a2) = 1.0;
        mapMatrix(3, 2 * a2 + 1) = 1.0;

        m_assembleStiffness += mapMatrix.transpose() *
                elem.m_K * mapMatrix;
    }

    m_globalStiffness = m_assembleStiffness;
    int j = 0;
    for(int index : indicesToConstraint){
        for(int row = 0; row < 2 * n; ++row){
            for(int col = 0; col < 2 * n; ++col){
                if (row == index || col == index){
                    m_globalStiffness(row, col) = 
                        (row == col) ? 1.0 : 0.0;
                }
            }
        }

        //removeRow(m_globalStiffness,index - j);
        //removeColumn(m_globalStiffness, index - j);
        j++;
    }
}

void trussStructure::recalcMass()
{
    int n = globalNodeList.size();
    m_assembleMass = Eigen::MatrixXd::Zero(2 * n, 2 * n);
    for(auto elem : trussArray){
        Eigen::MatrixXd mapMatrix = Eigen::MatrixXd::Zero(4, 2 * n);
        int a1 = elem.getNumNode1();
        int a2 = elem.getNumNode2();
        mapMatrix(0, 2 * a1) = 1.0;
        mapMatrix(1, 2 * a1 + 1) = 1.0;
        mapMatrix(2, 2 * a2) = 1.0;
        mapMatrix(3, 2 * a2 + 1) = 1.0;

        m_assembleMass += mapMatrix.transpose() *
                elem.m_mass * mapMatrix;

        std::cout << mapMatrix.transpose() *
                     elem.m_mass * mapMatrix << "\n\n";
    }
    m_globalMass = m_assembleMass;
    int j = 0;
    for(int index : indicesToConstraint){
        for(int row = 0; row < 2 * n; ++row){
            for(int col = 0; col < 2 * n; ++col){
                if (row == index || col == index){
                    m_globalMass(row, col) = 
                        (row == col) ? 1.0 : 0.0;
                }
            }
        }

        //removeRow(m_globalMass, index - j);
        //removeColumn(m_globalMass, index - j);
        j++;
    }
}

void trussStructure::recalcLoad()
{
    int n = globalNodeList.size();
    m_assembleLoad = Eigen::MatrixXd::Zero(2 * n, 1);
    for(int i = 0; i < n; ++i){
        m_assembleLoad(2 * i) += xLoadOnNode[i];
        m_assembleLoad(2 * i + 1) += yLoadOnNode[i];
    }

    m_globalLoad = m_assembleLoad;
    int j = 0;
    for(int index : indicesToConstraint){
        for(int row = 0; row < 2 * n; ++row){
                if (row == index){
                    m_globalLoad(row, 0) = 0.0;
                }
        }

        //removeRow(m_globalLoad, index - j);
        j++;
    }
}

这里留出比较方便的接口，例如建立一个结构的代码如下:

double l = .01;
double r = .005;
double a = 3.141592693975358 * r * r;
double e = 7e3 * 1e7;
double EA = e * a;
auto rho = 2778.62;

trussStructure structure1;
globalNodeList.push_back(node(0, 0));
globalNodeList.push_back(node(.5 * l, .5 * sqrt(3.0) * l));
globalNodeList.push_back(node(l, 0));

structure1.addTrussElement(trussElement(globalNodeList,
                                        0, 1, EA, rho));

structure1.addTrussElement(trussElement(globalNodeList,
                                        0, 2, EA, rho));

structure1.addTrussElement(trussElement(globalNodeList,
                                            1, 2, EA, rho));

structure1.addLoad(1, .0, 1e-3);
structure1.addConstraint(Constraint(1, 1, 0));
structure1.addConstraint(Constraint(0, 1, 2));


structure1.recalcAllMatrices();

就可以方便地建立一个桁架结构。不过还是需要进一步改善接口，增强其可扩展性。