MatLib 1.0 alpha

最新推荐文章于 2024-09-05 16:20:57 发布

原创最新推荐文章于 2024-09-05 16:20:57 发布 · 1.2k 阅读

0 ·

CC 4.0 BY-SA版权

文章标签：

#float #include

本文介绍了一种基于四节点单元的有限元分析方法，详细展示了如何通过计算局部刚度矩阵并将其组合成全局刚度矩阵来解决具体的力学问题。此外，还探讨了如何应用边界条件，并提供了一个具体的求解过程。

#include "../../MatLib/MathLib.h" #include <iostream> using namespace std; using namespace MathLib; using namespace MathLib::Matrix; int solve() { //和求解问题相关的常数 const Int32 nPoint = 12; const Int32 nDOF = 2; const Int32 nTotalDOF = nPoint * nDOF; const Int32 nFixedPoint = 2; const Int32 nElement = 5; const Int32 nNode = 4; const Int32 nEDOF = nNode * nDOF; // const Int32 nDime = 2; //材料常数 const Float young = 1500; const Float poiss = 0.01; //高斯积分常数 const Int32 nGauss = 2; const FullMatrix GaussWeight(2, 1, 1); FullMatrix GaussPoint(2, 1); GaussPoint.assigner() = -0.577350269189626, 0.577350269189626; //局部坐标到全局坐标的映射 FullMatrix LocalNodes(nElement, nNode, 0); { for (int part = 1; part <= nElement; ++part) { Int32 last = (part - 1) * 2; LocalNodes(part, 1) = last + 1; LocalNodes(part, 2) = last + 3; LocalNodes(part, 3) = last + 4; LocalNodes(part, 4) = last + 2; } } //全局坐标 FullMatrix Coord(12, 2, 0); { const Int32 Height = 1; Coord.assigner() = 0, 0, 0, Height, 1, 0, 2, Height, 2, 0, 4, Height, 4, 0, 5, Height, 7, 0, 6, Height, 10, 0, 10, Height; } //本质边界条件 const FullMatrix FixedValue(nTotalDOF, 1, 0); BasicMatrix<Int32> FixedPoint(1, nFixedPoint, 0); BasicMatrix<Int32> BoundDOF(nFixedPoint, nDOF, 0); BasicMatrix<Int32> IsFixed(1, nTotalDOF); { FixedPoint(1, 1) = 1; BoundDOF(1, 1) = 0; BoundDOF(1, 2) = 0; IsFixed(1, 1) = 1; IsFixed(1, 2) = 1; FixedPoint(1, 2) = 2; BoundDOF(2, 1) = 0; BoundDOF(2, 2) = 1; IsFixed(1, 3) = 1; IsFixed(1, 4) = 0; } //D矩阵 const Float E0 = young / (1 - poiss * poiss); const Float mu0 = poiss / (1 - poiss); FullMatrix DMatrix(3, 3, 0); { DMatrix.assigner() = 1, mu0, 0, mu0, 1, 0, 0, 0, 0.5 * (1 - mu0); DMatrix *= E0 / (1 - mu0 * mu0); } //总刚矩阵和总体载荷 FullMatrix GStiff(nTotalDOF, nTotalDOF, 0); FullMatrix GLoad(nTotalDOF, 1, 0); for (Int32 CurrElem = 1; CurrElem <= nElement; ++CurrElem) { //单元形状参数 FullMatrix XNode(nNode); FullMatrix YNode(nNode); for (Int32 i = 1; i <= nNode; ++i) { Int32 idx = (Int32)LocalNodes(CurrElem, i); XNode(i) = Coord(idx, 1); YNode(i) = Coord(idx, 2); } Float a1 = (-XNode(1) + XNode(2) + XNode(3) - XNode(4)) / 4; Float a2 = ( XNode(1) - XNode(2) + XNode(3) - XNode(4)) / 4; Float a3 = (-XNode(1) - XNode(2) + XNode(3) + XNode(4)) / 4; Float b1 = (-YNode(1) + YNode(2) + YNode(3) - YNode(4)) / 4; Float b2 = ( YNode(1) - YNode(2) + YNode(3) - YNode(4)) / 4; Float b3 = (-YNode(1) - YNode(2) + YNode(3) + YNode(4)) / 4; Float J0 = a1 * b3 - a3 * b1; Float J1 = a1 * b2 - a2 * b1; Float J2 = a2 * b3 - a3 * b2; //单刚矩阵 FullMatrix EStiff(nEDOF, nEDOF, 0); for (Int32 iGauss = 1; iGauss <= nGauss; ++iGauss) { Float x = GaussPoint(iGauss); for (Int32 jGauss = 1; jGauss <= nGauss; ++jGauss) { Float y = GaussPoint(jGauss); FullMatrix Deriv4(2, 4); { FullMatrix Right(4, 4); FullMatrix Left(1, 4); Left.assigner() = 0, 1, 0, y; Right.assigner() = 1, 1, 1, 1, -1, 1, 1,-1, -1,-1, 1, 1, 1,-1, 1,-1; Deriv4(1, 1, 1, 4) = 0.25 * Left * Right; Left.assigner() = 0, 0, 1, x; Deriv4(2, 2, 1, 4) = 0.25 * Left * Right; } Float djacb = J0 + J1 * x + J2 * y; Float dvolu = djacb * GaussWeight(iGauss) * GaussWeight(jGauss); FullMatrix xjaci(2, 2); { xjaci.assigner() = b3 + b2 * x, -(b1 + b2 * y), -(a3 + a2 * x), a1 + a2 * y; xjaci /= djacb; } FullMatrix bmatrix(3, 8); { FullMatrix cartd(2, 4); cartd = xjaci * Deriv4; bmatrix.assigner() = cartd(1, 1), 0, cartd(1, 2), 0, cartd(1, 3), 0, cartd(1, 4), 0, 0, cartd(2, 1), 0, cartd(2, 2), 0, cartd(2, 3), 0, cartd(2, 4), cartd(2, 1), cartd(1, 1), cartd(2, 2), cartd(1, 2), cartd(2, 3), cartd(1, 3), cartd(2, 4), cartd(1, 4); } EStiff += dvolu * Tran(bmatrix) * DMatrix * bmatrix; } } //单元载荷 FullMatrix EForce(nEDOF, 1, 0); { if(CurrElem == nElement) { EForce.assigner() = 0, 0, 1000, 0, -1000, 0, 0, 0; } } //单刚矩阵到总刚矩阵的组装 BasicMatrix<Int32> index(8); for (Int32 inode = 1; inode <= nNode; ++inode) { Int32 idx = (Int32)LocalNodes(CurrElem, inode); for (Int32 idof = 1; idof <= nDOF; ++idof) { index(nDOF * (inode - 1) + idof) = nDOF * (idx - 1) + idof; } } for (Int32 i = 1; i <= nEDOF; ++i) { Int32 idxi = index(i); GLoad(idxi) = GLoad(idxi) + EForce(i); for (Int32 j = 1; j <= nEDOF; ++j) { Int32 idxj = index(j); GStiff(idxi, idxj) += EStiff(i, j); } } } //使用本质边界条件 for (Int32 i = 1; i <= nFixedPoint; ++i) for (Int32 j = 1; j <= nDOF; ++j) { Int32 idx = (FixedPoint(1, i) - 1) * nDOF + j; if (IsFixed(1, idx)) { Float value = FixedValue(idx); for (Int32 itotv = 1; itotv <= nTotalDOF; ++itotv) { GLoad(itotv) -= value * GStiff(itotv, idx); GStiff(itotv, idx) = 0; GStiff(idx, itotv) = 0; } GStiff(idx, idx) = 1; GLoad(idx, 1) = value; } } //求解 FullMatrix res(GStiff.row(), 1); GaussSolve(res, GStiff, GLoad); cout.precision(16); cout<< res; cout << res(21) << '/t' << res(22) << endl; return 0; }

给我发邮件可以获得本库的代码

mail:bailiang@vip.qq.com