等几何分析+有限元等参四边形单元的等效节点力

---

前言

写等几何程序的时候用到等效节点力，只推导了2次B样条单元（一个单元九个节点）的等效节点力公式。

一、有限元等参四边形单元的等效节点力

1.1 有限元等参四边形单元理论

可参考博客：等参单元与数值积分 或者图书《有限元法与MATLAB程序设计》

1.2 等参四边形单元的等效节点力

内容来自于：图书《有限元法与MATLAB程序设计》

二、等几何分析的等效节点力

2.1 等几何分析的理论

等几何分析的相关理论可以查看博客：等几何分析编程记录

2.2 等效节点力

理论推导：

写程序时的分析：

对应的C++函数：

// 求等效节点载荷公式  在参数域  的函数 N' *  P * |jacobi|
	MatrixXd return_F_function(MatrixXd WF, MatrixXd elCpts, double u, int num_of_ele)
	{
		// ------------------------------------------------------------------------------------------------------
		// 计算边界的等效节点载荷公式中积分里面的函数 N' *  P * |jacobi|
		// 输入
		//      u:变换到标准高斯区间的高斯点
		//     WF:外力载荷
		//         2维:  WF = [Px, Py] 
		//         3维:  WF = [Px, Py, Pz] 
		//     elCpts : 控制点坐标 ------- 将边界点坐标传入 从下到上（对应基函数）
		//              coordinates of element control points 单元控制点坐标
		// 输出: 
		//   F --- 等效节点载荷
		//         不细化  F:6行1列的向量
		// -------------------------------------------------------------------------------------------------------
		//  边界对应的形函数（基函数）
		MatrixXd N(2, 6);
		N.setZero();
		int len_V = size(V);
		int num1 = num_of_ele;
		for (int i = 0; i < 3; i++)
		{
			N(0, 2 * i) = OneBasisFun_withreturn(2, len_V - 1, V, num1, u);
			N(1, 2 * i + 1) = OneBasisFun_withreturn(2, len_V - 1, V, num1, u);
			num1++;
		}

		//  边界对应的形函数导数（基函数）
		MatrixXd der_B(1, 3);
		//Print_Vector(der_B);
		int num2 = num_of_ele;
		for (int i = 0; i < 3; i++)
		{
			der_B(0, i) = DersOneBasisFun_double_First_Ders(2, V, num2, u);
			num2++;
		}
		double j2 = 0; // jacobi_2：物理域到参数域的jacobi值
		j2 = (der_B * elCpts).norm(); // der_B * elCpts:导函数 * 对应的控制点坐标

		//MatrixXd F = t * j2 * N.transpose() * WF;
		MatrixXd F = j2 * N.transpose() * WF;
		return F;
	}

	// 求单元的等效节点载荷
	MatrixXd ele_F(vector<double> elDoma, MatrixXd WF, MatrixXd elCpts, vector<int> Gauss_num, int num_of_ele)
	{
		// ------------------------------------------------------------------------------------------------------
		// 计算边界单元的等效节点载荷
		// 目前适用于只有右边界施加力的等效节点载荷
		// 输入
		//     elDoma : 传入的单元参数域     
		//         1维:   elDoma = [u1, u2]     
		//         2维:   elDoma = [u1, u2, v1, v2] 
		//         3维:   elDoma = [u1, u2, v1, v2, w1, w2] 
		//     WF:应力
		//         2维:  WF = [Px, Py] 
		//         3维:  WF = [Px, Py, Pz] 
		//     elCpts : 控制点坐标 ------- 将边界点坐标传入 从下到上（对应基函数）
		//   Gauss_num: 高斯积分用到的高斯点数
		//              因为是一维高斯积分，故Gauss_num={2} ----  表示用两点高斯积分求解     
		//              或者Gauss_num={3} ----  表示用3点高斯积分求解 
		// 输出: 
		//   f --- 单元等效节点载荷 
		// -------------------------------------------------------------------------------------------------------
		Gauss gau; //定义一个高斯类，用于调用"gauss.h"中函数
		MatrixXd f(6, 1);
		f.setZero();

		gau.gauss_quadrature(Gauss_num);  //必须先调用这个函数，才能得到 gp 和 wgt 的值
		MatrixXd gp = gau.standard_GaussPoint; //标准高斯区间的高斯点
		int num_gp = gp.cols();//高斯积分用到的高斯点数
		MatrixXd wgt = gau.standard_GaussWeight;//标准高斯区间的权值
		//cout << " gp = \n" << gp << endl;  cout << " wgt = \n" << wgt << endl;
		double u = 0;
		double j1 = gau.jacobian_gauss_mapping(elDoma); //高斯积分映射---参数域区间到标准高斯区间的jacobi

		for (int i = 0; i < num_gp; i++)
		{
			double a = elDoma[0];  //积分下限
			double b = elDoma[1];  //积分上限
			u = (b - a) / 2 * gp(0, i) + (a + b) / 2;
			f += j1 * wgt(i) * return_F_function(WF, elCpts, u, num_of_ele);
		}
		//cout << " f = \n" << f << endl;
		return f;
	}

	//计算整体等效节点载荷
	MatrixXd return_F(VectorXd val_WF)
	{
		// ------------------------------------------------------------------------------------------------------
		// 计算整体等效节点载荷
		// 目前适用于只有右边界施加力的等效节点载荷
		// 输入;
		//       rightNodes: 右边界节点的全局编号
		//     num_of_ele_V: V方向单元个数                                       
		//        Gauss_num: 高斯积分用到的高斯点数   Gauss_num={2} ----  表示用两点高斯积分求解
		//           val_WF: 外力的大小    2维:  WF = [Px, Py]                            
		//          elCpts : 边界控制点坐标 ------- 将边界点坐标传入 从下到上（对应基函数）                                
		// 输出: 
		//             F --- 等效节点载荷        
		//                             
		//   注:  2次基函数, V方向一个单元控制点 只有三个                                          
		// -------------------------------------------------------------------------------------------------------

		MatrixXd F(num_GK, 1);
		F.setZero();
		vector<int> Gauss_num = { 4 }; //高斯积分用到的高斯点数
		//计算施加外力所对应节点的坐标
		num_BasisFun_V;  // V方向基函数个数
		MatrixXd Cpts(num_BasisFun_V, 2);   // elCpts:coordinates of control points 边界控制点坐标
		Cpts.setZero();
		//外力施加在右边界
		for (int i = 0; i < num_BasisFun_V; i++)
		{
			Cpts.row(i) = Node.row(rightNodes[i]);//右边界节点的全局编号对应的坐标取出来
		}
		//elCpts.block<rows, cols>(i, j)表示 从矩阵elCpts的  elCpts(i, j) 位置开始取rows行  cols列
		//计算单元等效节点力
		double Area = t * L; //截面积
		val_WF /= Area;//应力
		/*cout << " val_WF = \n" << val_WF << endl;
		cout << "\n 力边界节点: rightNodes = \n"; Print_Vector(rightNodes);*/

		//组装等效节点力
		for (int num = 0; num < num_of_ele_V; num++)
		{
			// elCpts:coordinates of element control points 单元控制点坐标
			MatrixXd elCpts = Cpts.block<3, 2>(num, 0); //取出右边界（V方向）的 第num个单元的坐标
			// 计算单元节点坐标  对应的参数域
			// elDoma:参数域的积分区间
			vector<double> elDoma = { val_V[num],val_V[num + 1] };
			MatrixXd f1 = ele_F(elDoma, val_WF, elCpts, Gauss_num, num);//求出第num个单元的等效节点力
			//cout << " f1 = \n" << f1 << endl;
			vector<int> num_ele = { rightNodes[num], rightNodes[num + 1],rightNodes[num + 2] };//第num个单元对应节点的全局编号
			//cout << "\n 力边界节点: num_ele = \n";Print_Vector(num_ele);

			//组装等效节点力载荷  将单元等效节点力放到相应位置
			//注意:节点从0开始编号
			//若节点num_ele[0] = 5  对应在等效节点力向量中 2*5 = 10 和 11 行的位置
			F(2 * num_ele[0]) += f1(0);
			F(2 * num_ele[0] + 1) += f1(1);
			F(2 * num_ele[1]) += f1(2);
			F(2 * num_ele[1] + 1) += f1(3);
			F(2 * num_ele[2]) += f1(4);
			F(2 * num_ele[2] + 1) += f1(5);
		}
		//cout << " F = \n" << F << endl;
		return F;
	}