询问一个关于高斯点应力计算的问题

shikang

       在有限元计算当中，单元内高斯点的应力精度要高于节点上的应力。也因此Abaqus一般给出的应力是高斯点上的应力。如果描述有误，请指正。下面将正式询问一个问题：之所以会存在高斯点，因为在采用等参元计算单元刚度矩阵时，采用了高斯积分。这里以一个简单的六面体单元为例，如果在积分时，每个方向的积分选择2个高斯点，那么在这个单元内就会有2×2×2=8个高斯点，如果每个方向选择3个高斯点，则这个单元内会有3×3×3=27个高斯点。不管选择多少个高斯点，对于同一个单元，计算好的单元刚度矩阵维数是一定的。那么问题就来了，如果我通过应力计算公式[S]=[D]×[B]×[d]得到的到底是节点上的应力还是高斯点的应力呢？如果是节点应力，那如何计算高斯点的应力呢（注意，同一个单元可能有8个高斯点也可能有27个高斯点）？如果是高斯点应力，但是如果我采用的27个高斯点计算，那这个维度似乎是不正确的       编程的路很孤独，希望得到您的指点？

penultimate

[S]是节点应力，不是高斯点。找个简单点儿的单元，推导一下单刚矩阵就清楚了

hillyuan

1。一般情况下没得不到节点应力，你需要从高斯点外插。
2。如果应力在三方向上都是线性分布，8个高斯点足够，采用27个高斯点的计算结果理论上与8个高斯点相同
3。数值计算精度问题，与维度无关

penultimate

penultimate 发表于 2015-7-28 20:27
[S]是节点应力，不是高斯点。找个简单点儿的单元，推导一下单刚矩阵就清楚了 ...

原始的[D]*[B]*[d]相乘的结果的确是高斯点处的应力，但对于最后的应力平滑，求高斯点处的应力没有意义，应该外延到节点，求节点处的应力。
谢谢hillyuan指出我的错误

shikang

hillyuan 发表于 2015-7-28 21:31
1。一般情况下没得不到节点应力，你需要从高斯点外插。
2。如果应力在三方向上都是线性分布，8个高斯点足够 ...

虽然应力在三个方向是线性分布，但这时我却采用了27个高斯点去计算，那么按照应力计算公式却只能计算出8个点应力的值（这时的应力值和采用8个高斯点计算出来的值差不多），而不是计算出27个，如果这8个点的应力是高斯点的应力，那这时的高斯点和27个高斯点有何区别呢？

骞耕

d是节点位移矩阵，是个常数阵，雷打不动；B是形函数偏导阵，形函数是自然坐标的函数，把哪点的自然坐标代入，[s]=[D][B][d]求的就是哪点的应力。这样回答不知道我有没有正确理解楼主的意思。

骞耕

把节点处的自然坐标代入公式 [s]=[D][B][d]  也是可以求节点应力的，但这样做没有先算高斯点再外插精确。

骞耕

对于线性单元，8个高斯点精度足够，当然也可以使用27个高斯点。把27个高斯点自然坐标代入 [s]=[D][B][d] ，可以求出27个高斯点的应力。如何根据27个高斯点应力外推得到8个节点应力？先说下一般的外推方法（见王勖成老师的书）：节点结果已知，使用插值函数可以得到场内任一点值， 高斯点的值也可以通过这种方法得到，即 [N]*[节点值]=[高斯点值]，但我们是先得到高斯点值的，好办[N逆]*[高斯点值]=[节点值]，如此便将高斯点值反推到了节点；[N]可求逆则[N]必为方阵，但对于8节点、27个高斯点的单元而言[N]的维度为27*8（行为高斯点，列为节点），不是方阵，怎么办？将单元的插值函数由8个扩展为27个，则[N]变成了27*27，通过 [N逆]*[高斯点值]=[节点值] 得到27个节点值，去掉增加的只保留8个就可以了。

shikang

骞耕 发表于 2015-7-29 10:18
把节点处的自然坐标代入公式 =[D][d]  也是可以求节点应力的，但这样做没有先算高斯点再外插精确。 ...

1、很感谢回答。
2、如果是这样，那分歧点就没有了。在这里说一下最后的理解。我主要是把前后两个B矩阵给混淆了。在组装刚度矩阵时，需要求解B矩阵，而这时的B矩阵是通过节点坐标求得的。在计算应力时=[D][B][d]，也有一个B矩阵，而这时的B矩阵代入不同的点的坐标就得到不同点节点的应力，代入节点坐标得到对应节点的应力，代入高斯点的坐标得到高斯点的应力，而高斯点坐标的应力值比较精确。
3、对于从高斯点应力推算节点应力，就是采用合适的插值算法求得。

感谢上面所有人的回答！

骞耕

shikang 发表于 2015-7-29 18:24
1、很感谢回答。
2、如果是这样，那分歧点就没有了。在这里说一下最后的理解。我主要是把前后两个B矩阵给 ...

其实。。。。。。
组装刚度矩阵时，计算的[B]也是高斯点上的。

QustFelix

骞耕 发表于 2015-7-30 15:39
其实。。。。。。
组装刚度矩阵时，计算的也是高斯点上的。

我也有点疑问：
1.单元内高斯点的应力精度要高于节点上的应力，我还没有读过《有限元方法》一书，可否简要说明一下为什么?
2.我理解一阶三角形平面应力单元，[B]是个常数矩阵，应变【ε】 = [B]*[δ], B为3*6 矩阵，δ为6×1，所以ε为3×1，所以一阶平面单元为常应变，也就是常应力单元。
但是对于二阶三角形平面应力单元，[B]是单元坐标的函数，应变【ε】 = [B]*[δ], B为3*12 矩阵，δ为12×1，所以ε为3×1，但是对于不同点，有不同的B，所以有不同的ε，也就是单元内各点应力也不同。也就意味着二阶单元是非常应力的。比方说我取单元内的三个高斯点，当然可以计算出这三个高斯点的应变，应力。

我的理解是否正确呢？ 谢谢指点！

骞耕

关于您的疑问，做以下回答：
      1.我也不知道为什么高斯点精确，书上好像也没有解释，大概是因为组装单刚时采用高斯积分的缘故；
      2.您的理解和正确。

快乐宙斯

如果应力在三方向上都是线性分布，8个高斯点足够，采用27个高斯点的计算结果理论上与8个高斯点相同

refeihc

同意6、7楼的观点，B就是几何矩阵，一般情况下随坐标变化，代入高斯点的坐标就得到该点的相应值。

至于高斯点的应力精度要高于节点上的应力，这是基于前人超收敛性研究的成果总结出来的，我等一般可抛开严格的证明过程不去管，但是要注意几个要点。

1  若场函数采用n次多项式逼近，其精度可达到n+1；
2  应力和应变是场函数求一阶导数，为n-1次多项式，所以单元上应力和应变的精度就要降为n阶；
3  应力和应变的近似解会在真实解附近振荡，因此近似曲线（或曲面、超曲面）和精确曲线必然在一些地方有相交，虽然理论上不能准确求出，但高斯积分点很靠近这些相交的地方，这也就是在高斯积分点处应力和应变的精度会高一阶的原因。

refeihc

针对3、5、13楼的讨论，谈点看法。

所谓应力在三个方向是线性分布，这一说法很奇怪，有限元数值计算中有什么理由单把应力做为被积函数呢？再谈及因而采用8个和27个高斯积分点的效果相同，这一说法没有意义。
猜想大家本来是想说位移插值函数吧，或者是单元刚度矩阵和单元节点荷载向量的数值计算。

一般情况下，8结点六面体等参元的位移插值函数对局部坐标的三个方向的确是呈线性分布，但当单元形状不规则时，即使是线弹性问题，被积函数也未必在三个局部坐标方向上呈线性分布。

只要想一想应变矩阵的计算中涉及到对全局坐标求偏导，会出现逆矩阵，以及求面积分时，面积元素因坐标变换会出现根式，就能理解了。

所以采用8个高斯点和27个高斯点得到的计算结果一般是不相同的。

只采用8个高斯积分点的原因是因为精度够了，虽然更多积分点会带来数值积分精度的提高，但是于整体精度的提高无补，反而造成过大的计算开销。

最后补充说一句，对于六面体单元，采用个高斯积分点获得的积分效果未必是最优的，有更节省计算成本的方案，可能是因为高斯积分太受编程者的欢迎了吧，在数值积分这个环节上再进一步提高效率的动机变弱了。