高斯积分相关

​ 对于一般的积分表达式，是无法求出其定积分形式的。对于有限元问题而言，主要就是单元刚度矩阵的问题。（当然对于梁单元而言，是可以给出单元刚度矩阵显式表达式的）。因此得用到数值积分方式。而这其中，如果积分具有如下形式，则可以应用高斯积分。

$$\int^1_{-1}f(x)\text{d}x\\ \int^1_{-1}\int^1_{-1}f(x,y)\text{d}x\text{d}y\\ \int^1_{-1}\int^1_{-1}\int^1_{-1}f(x,y,z)\text{d}x\text{d}y\text{d}y\\ $$

高斯积分表达式为： $$\sum^n_{i=1}H_if(x_i)\\ \sum^n_{i=1}\sum^n_{i=1}H_{ij}f(x_{ij},y_{ij})\\ \sum^n_{k=1}\sum^n_{j=1}\sum^n_{i=1}H_{ijk}f(x_{ijk},y_{ijk},z_{ijk})\\ $$ 可以看出，高斯积分实质上可视为一种加权点积分的形式，其优点在于如果被积函数是多项式的话，高斯积分是可以做到准确的，即有：

对于n个积分点，只要选取适当的加权系数和积分点位置，能够使积分式在被积函数为不超过（2n-1）次多项式时精确成立。

可以看出实际应用高斯积分时，需要满足两个条件：

（1）即各维度积分上下限为1，-1，这个可以通过参数变换实现。

（2）被积函数为多项式

​ 对于平面问题形心 $\int{x}\text{d}A$，惯性矩 $\int{x^2}\text{d}A$，对于满足平截面假定的弯矩求解等均可以适用。

对于一般的有限元而言，单元刚度矩阵形式可表达为： $$\int_V{B^\rm T DB}\text{d}V$$ 经等参单元表达后积分形式可统一表示为： $$\int^1_{-1}\int^1_{-1}\int^1_{-1}f(\xi,\eta,\gamma)\text{d}\xi\text{d}\eta\text{d}\gamma=\int^1_{-1}\int^1_{-1}\int^1_{-1}{B(\xi,\eta,\gamma)^\rm T D(\xi,\eta,\gamma)B(\xi,\eta,\gamma)}|J|\text{d}\xi\text{d}\eta\text{d}\gamma$$ ​ 对于有限元计算而言，常规的梁、杆单元而言，由于刚度矩阵具有明确的解析表达式，因而也就不会有数值积分这一问题了。但对于而对于目前的有限元方法主要是采用等参单元以及非协调元，由于子单元非常复杂，通常情况下导致Jacobi及其行列式计算比较复杂，此时一般都不能进行显示积分计算(刚度矩阵$K$的计算中含有雅可比行列式)，而需要借助于数值积分方法。而能否采用高斯积分关键在于$f(\xi,\eta,\gamma)$是否为多项式。

​ 不失一般性，以通常的四边形单元来说明，采用拉格朗日插值可以表示为: $$x=\left(\frac{\xi-1}{-2}\right)\left(\frac{\eta-1}{-2}\right)x_{11}+\left(\frac{\xi+1}{2}\right)\left(\frac{\eta-1}{-2}\right)x_{12}+\left(\frac{\xi-1}{-2}\right)\left(\frac{\eta+1}{2}\right)x_{21}+\left(\frac{\xi+1}{2}\right)\left(\frac{\eta+1}{2}\right)x_{22}\\ y=\left(\frac{\xi-1}{-2}\right)\left(\frac{\eta-1}{-2}\right)y_{11}+\left(\frac{\xi+1}{2}\right)\left(\frac{\eta-1}{-2}\right)y_{12}+\left(\frac{\xi-1}{-2}\right)\left(\frac{\eta+1}{2}\right)y_{21}+\left(\frac{\xi+1}{2}\right)\left(\frac{\eta+1}{2}\right)y_{22}$$ 基于上述形式，雅可比矩阵$|J|$为关于$\xi,\eta$的二次多项式。又形函数$N$通常为多项式形式，故$B^\rm T DB$也一般为多项式，用$\xi,\eta$变量代换后亦为多项式，故

$f(\xi,\eta)$应为关于$\xi,\eta$的多项式，即可以采用高斯积分。在这篇ppthttps://wenku.baidu.com/view/1ce457b6a32d7375a5178043.html中对于有限元中高斯积分应用有着详细的介绍。

​ 若单元为矩形， 则： $$x=\xi (w/2)cos\theta-\eta (h/2)sin\theta+x_c\\ y=\xi (w/2)sin\theta+\eta (h/2)cos\theta+y_c$$ 这里$w、h$分别为矩形宽度和高度，$x_c、y_c$分别为矩形中心点在整体坐标系中的x、y坐标,$\theta$w为矩形转角。相应的，雅可比行列式为： $$\left|J\right|=\left|\begin{array}{}(w/2)cos\theta&-(h/2)sin\theta\\(w/2)sin\theta&(h/2)cos\theta\end{array}\right|=wh/4$$ ​ 然后有个缩减积分概念,即对于2n-1次多项式，如果采用的积分点少于n，那么理论上会无法使积分完全一致，但可以提高计算效率。

​ 高斯积分的另一大问题在于计算结果的外推（也即插值），上述这个ppt中提到：

应力一般采用多个积分点的相互插值或外延来计算节点应力。这只是为了减少误差。因为在积分点应力比节点具有更高阶的误差。

​ 对于Beam188单元，keyopt(3)用于设置形函数类型。如果keyopt(3)=0默认情形，只有一个积分点，相当于在外推时只有一个插值点，对于单元的量值都是均一的。那么所有单元结果量将在全长范围内保持一致。换句话说，弯矩将是保持一致。只有纯弯情形下才是对的。故一般需要进行更为精细的单元划分。经过验证，剪力不是基于弯矩导数的，而是单独计算的，故也沿全长相等。如果keyopt(3)=2，那么在将有一内部节点。共有两个高斯积分点。相应所有单元结果量将在全长范围内呈现线性。

相关程序:

# 二维9积分点高斯积分，可适用于最高次数为5次的平面积分
import numpy as np

# 各积分点在局部坐标系中的坐标位置
dg= 0.774596669241483
Guass9_X = np.array([[-dg,-dg,-dg],[0,dg,dg],[dg,0,0]])
Guass9_Y = np.array([[-dg,0,dg],[dg,dg,0],[-dg,-dg,0]])

# 各积分点的权重
we1 = 0.308641975308642
we2 = 0.493827160493828
we3 = 0.790123456790124
weight9 = np.array([[we1,we2,we1],[we2,we1,we2],[we1,we2,we3]])

def Gauss9f(f):   
    ret = np.sum(weight9*f)
    return ret 

对于二维高斯积分点及其权重，可参见http://blog.sina.com.cn/s/blog_625936040100qs3d.html。

2020年8月24日 星期一 天气晴