20 结点四面体单元的形函数

写出图示20结点四面体单元的形函数（选做）解：记四面体单元的体积为V，A(x,y,z)为单元内任一点，1、2、3、4的坐标分别为(x1,y1,z1)、(x2,y2,z2)、(x3,y3,z3)、(x4,y4,z4)为A点的体积坐标由本次作业第1题知，,其中a,b,c,d参数同课本P60下面利用几何法构造形函数：对结点1，作3个平面过除1外的所有节点但不过1，即：则对结点2，作3个平面过除2外的所有节点但不过2，即：则依此类推可得到：N3=1/2*L3*(3*L3-1)*(3*L3-2);N4=1/2*L4*(3*L4-1)*(3*L4-2);N5=9/2*L1*L2*(3*L1-1);N6=9/2*L1*L2*(3*L2-1);N7=9/2*L1*L4*(3*L1-1);N8=9/2*L1*L4*(3*L4-1);N9=9/2*L1*L3*(3*L1-1);N10=9/2*L1*L3*(3*L3-1);N11=9/2*L4*L3*(3*L3-1);N12=9/2*L4*L3*(3*L4-1);N13=9/2*L4*L2*(3*L4-1);N14=9/2*L4*L2*(3*L2-1);N15=9/2*L3*L2*(3*L3-1);N16=9/2*L3*L2*(3*L2-1);N17=27*L1*L2*L4;N18=27*L3*L2*L4;N19=27*L1*L2*L3;N20=27*L1*L3*L4;接下来进行完备性和协调性检验，利用matlab编写如下程序：symsxyzx1x2x3x4y1y2y3y4z1z2z3z4%结点坐标x5=2/3*x11/3*x2;y5=2/3*y11/3*y2;z5=2/3*z11/3*z2;x6=2/3*x21/3*x1;y6=2/3*y21/3*y1;z6=2/3*z21/3*z1;x7=2/3*x11/3*x4;y7=2/3*y11/3*y4;z7=2/3*z11/3*z4;x8=2/3*x41/3*x1;y8=2/3*y41/3*y1;z8=2/3*z41/3*z1;x9=2/3*x11/3*x3;y9=2/3*y11/3*y3;z9=2/3*z11/3*z3;x10=2/3*x31/3*x1;y10=2/3*y31/3*y1;z10=2/3*z31/3*z1;x11=2/3*x31/3*x4;y11=2/3*y31/3*y4;z11=2/3*z31/3*z4;x12=2/3*x41/3*x3;y12=2/3*y41/3*y3;z12=2/3*z41/3*z3;x13=2/3*x41/3*x2;y13=2/3*y41/3*y2;z13=2/3*z41/3*z2;x14=1/3*x42/3*x2;y14=1/3*y42/3*y2;z14=1/3*z42/3*z2;x15=2/3*x31/3*x2;y15=2/3*y31/3*y2;z15=2/3*z31/3*z2;x16=1/3*x32/3*x2;y16=1/3*y32/3*y2;z16=1/3*z32/3*z2;x17=1/3*(x1x2x4);y17=1/3*(y1y2y4);z17=1/3*(z1z2z4);x18=1/3*(x3x2x4);y18=1/3*(y3y2y4);z18=1/3*(z3z2z4);x19=1/3*(x1x2x3);y19=1/3*(y1y2y3);z19=1/3*(z1z2z3);x20=1/3*(x1x3x4);y20=1/3*(y1y3y4);z20=1/3*(z1z3z4);%体积坐标V=1/6*det([1,x1,y1,z1;1,x2,y2,z2;1,x3,y3,z3;1,x4,y4,z4]);L1=1/(6*V)*det([1,x,y,z;1,x2,y2,z2;1,x3,y3,z3;1,x4,y4,z4]);L2=-1/(6*V)*det([1,x,y,z;1,x1,y1,z1;1,x3,y3,z3;1,x4,y4,z4]);L3=1/(6*V)*det([1,x,y,z;1,x1,y1,z1;1,x3,y3,z3;1,x2,y2,z2]);L4=1-L1-L2-L3;%L4=-1/(6*V)*det([1,x,y,z;1,x1,y1,z1;1,x2,y2,z2;1,x3,y3,z3]);%形函数N1=1/2*L1*(3*L1-1)*(3*L1-2);N2=1/2*L2*(3*L2-1)*(3*L2-2);N3=1/2*L3*(3*L3-1)*(3*L3-2);N4=1/2*L4*(3*L4-1)*(3*L4-2);N5=9/2*L1*L2*(3*L1-1);N6=9/2*L1*L2*(3*L2-1);N7=9/2*L1*L4*(3*L1-1);N8=9/2*L1*L4*(3*L4-1);N9=9/2*L1*L3*(3*L1-1);N10=9/2*L1*L3*(3*L3-1);N11=9/2*L4*L3*(3*L3-1);N12=9/2*L4*L3*(3*L4-1);N13=9/2*L4*L2*(3*L4-1);N14=9/2*L4*L2*(3*L2-1);N15=9/2*L3*L2*(3*L3-1);N16=9/2*L3*L2*(3*L2-1);N17=27*L1*L2*L4;N18=27*L3*L2*L4;N19=27*L1*L2*L3;N20=27*L1*L3*L4;x0=[x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20];y0=[y1,y2,y3,y4,y5,y6,y7,y8,y9,y10,y11,y12,y13,y14,y15,y16,y17,y18,y19,y20];z0=[z1,z2,z3,z4,z5,z6,z7,z8,z9,z10,z11,z12,z13,z14,z15,z16,z17,z18,z19,z20];N0=[N1,N2,N3,N4,N5,N6,N7,N8,N9,N10,N11,N12,N13,N14,N15,N16,N17,N18,N19,N20];%完备性和协调性检验A=simplify(N1N2N3N4N5N6N7N8N9N10N11N12N13N14N15N16N17N18N19N20)B=simplify(x0*N0')C=simplify(y0*N0')D=simplify(z0*N0')最终输出结果为：>>A=1B=1C=1D=1此即是说明上述形函数形式满足完备性和协调性，即它们是20结点四面体单元的形函数综上可得，20结点四面体单元的形函数为：N1=1/2*L1*(3*L1-1)*(3*L1-2);N2=1/2*L2*(3*L2-1)*(3*L2-2);N3=1/2*L3*(3*L3-1)*(3*L3-2);N4=1/2*L4*(3*L4-1)*(3*L4-2);N5=9/2*L1*L2*(3*L1-1);N6=9/2*L1*L2*(3*L2-1);N7=9/2*L1*L4*(3*L1-1);N8=9/2*L1*L4*(3*L4-1);N9=9/2*L1*L3*(3*L1-1);N10=9/2*L1*L3*(3*L3-1);N11=9/2*L4*L3*(3*L3-1);N12=9/2*L4*L3*(3*L4-1);N13=9/2*L4*L2*(3*L4-1);N14=9/2*L4*L2*(3*L2-1);N15=9/2*L3*L2*(3*L3-1);N16=9/2*L3*L2*(3*L2-1);N17=27*L1*L2*L4;N18=27*L3*L2*L4;N19=27*L1*L2*L3;N20=27*L1*L3*L4;其中L1，L2，L3，L4为体积坐标