有限元编程的c++实现算例
有限元编程的c++实现算例read.pudn./downloads76/doc/fileformat/290377/ganjian.cpp__.htm1. #include<stdio.h>2. #include<math.h>3.4.5. #define ne 3//单元数6. #define nj 4//节点数7. #define nz 6//支撑数8. #define npj 0//节点载荷数9. #define npf 1//非节点载荷数10. #define nj3 12//节点位移总数11. #define dd 6//半带宽12. #define e0 2.1E8//弹性模量13. #define a0 0.008//截面积14. #define i0 1.22E-4//单元惯性距15. #define pi 3.14159265416.17.18. int jm[ne+1][3]={{0,0,0},{0,1,2},{0,2,3},{0,4,3}};/*gghjghg*/19. double gc[ne+1]={0.0,1.0,2.0,1.0};20. double gj[ne+1]={0.0,90.0,0.0,90.0};21. double mj[ne+1]={0.0,a0,a0,a0};22. double gx[ne+1]={0.0,i0,i0,i0};23. int zc[nz+1]={0,1,2,3,10,11,12};24. double pj[npj+1][3]={{0.0,0.0,0.0}};25. double pf[npf+1][5]={{0,0,0,0,0},{0,-20,1.0,2.0,2.0}};26. double kz[nj3+1][dd+1],p[nj3+1];27. double pe[7],f[7],f0[7],t[7][7];28. double ke[7][7],kd[7][7];29.30.31. //**kz[][]—整体刚度矩阵32. //**ke[][]—整体坐标下的单元刚度矩阵33. //**kd[][]—局部坐标下的单位刚度矩阵34. //**t[][]—坐标变换矩阵35.36. //**这是函数声明37. void jdugd(int);38. void zb(int);39. void gdnl(int);40. void dugd(int);41.42.43. //**主程序开始44. void main()45. {46.int i,j,k,e,dh,h,ii,jj,hz,al,bl,m,l,dl,zl,z,j0;47.double cl,wy[7];48.int im,in,jn;49.50. //***********************************************51. //<功能：形成矩阵P>52. //***********************************************53.54.if(npj>0)55.{56.for(i=1;i<=npj;i++)57.{//把节点载荷送入P58.j=pj[i][2];59.p[j]=pj[i][1];60.}61.}62.if(npf>0)63.{64.for(i=1;i<=npf;i++)65.{//求固端反力F066.hz=i;67.gdnl(hz);68.e=(int)pf[hz][3];69.zb(e);//求单元70.for(j=1;j<=6;j++)//求坐标变换矩阵T71.{72.pe[j]=0.0;73.for(k=1;k<=6;k++)//求等效节点载荷74.{75.pe[j]=pe[j]-t[k][j]*f0[k];76.}77.}78.al=jm[e][1];79.bl=jm[e][2];80.p[3*al-2]=p[3*al-2]+pe[1];//将等效节点载荷送到P中81.p[3*al-1]=p[3*al-1]+pe[2];82.p[3*al]=p[3*al]+pe[3];83.p[3*bl-2]=p[3*bl-2]+pe[4];84.p[3*bl-1]=p[3*bl-1]+pe[5];85.p[3*bl]=p[3*bl]+pe[6];86.}87.}88.89.90.//*************************************************91.//<功能：生成整体刚度矩阵kz[][]>92.for(e=1;e<=ne;e++)//按单元循环93.{94.dugd(e);//求整体坐标系中的单元刚度矩阵ke95.for(i=1;i<=2;i++)//对行码循环96.{97.for(ii=1;ii<=3;ii++)98.{99.h=3*(i-1)+ii;//元素在ke中的行码100.dh=3*(jm[e][i]-1)+ii;//该元素在KZ中的行码101.for(j=1;j<=2;j++)102.{103.for(jj=1;jj<=3;jj++)//对列码循环104.{105.l=3*(j-1)+jj;//元素在ke中的列码106.zl=3*(jm[e][j]-1)+jj;//该元素在KZ中的行码107.dl=zl-dh+1;//该元素在KZ*中的行码108.if(dl>0)109.kz[dh][dl]=kz[dh][dl]+ke[h][l];//刚度集成110.}111.}112.}113.}114.}115.116. //**引入边界条件**117.for(i=1;i<=nz;i++)//按支撑循环118.{119.z=zc[i];//支撑对应的位移数120.kz[z][l]=1.0;//第一列置1121.for(j=2;j<=dd;j++)122.{123.kz[z][j]=0.0;//行置0124.}125.if((z!=1))127.if(z>dd)128.j0=dd;129.else if(z<=dd)130.j0=z;//列（45度斜线）置0 131.for(j=2;j<=j0;j++)132.kz[z-j+1][j]=0.0;133.}134.p[z]=0.0;//P置0135.}136.137.138.139.140.for(k=1;k<=nj3-1;k++)141.{142.if(nj3>k+dd-1)//求最大行码143.im=k+dd-1;144.else if(nj3<=k+dd-1)145.im=nj3;146.in=k+1;147.for(i=in;i<=im;i++)148.{149.l=i-k+1;150.cl=kz[k][l]/kz[k][1];//修改KZ151.jn=dd-l+1;152.for(j=1;j<=jn;j++)153.{154.m=j+i-k;155.kz[i][j]=kz[i][j]-cl*kz[k][m];156.}157.p[i]=p[i]-cl*p[k];//修改P159.}160.161.162.163.164.p[nj3]=p[nj3]/kz[nj3][1];//求最后一个位移分量165.for(i=nj3-1;i>=1;i--)166.{167.if(dd>nj3-i+1)168.j0=nj3-i+1;169.else j0=dd;//求最大列码j0170.for(j=2;j<=j0;j++)171.{172.h=j+i-1;173.p[i]=p[i]-kz[i][j]*p[h];174.}175.p[i]=p[i]/kz[i][1];//求其他位移分量176.}177.printf("\n");178.printf("_____________________________________________________________\n");179.printf("NJ U V CETA\n");//输出位移180.for(i=1;i<=nj;i++)181.{182.printf(" %-5d %14.11f%14.11f%14.11f\n",i,p[3*i-2],p[3*i-1],p[3*i]); 183.}184.printf("_____________________________________________________________\n");185.//*根据E的值输出相应E单元的N,Q,M(A,B)的结果**186.printf("E N Q M\n");187.//*计算轴力N,剪力Q,弯矩M*188.for(e=1;e<=ne;e++)//按单元循环189.{190.jdugd(e);//求局部单元刚度矩阵kd191.zb(e);//求坐标变换矩阵T192.for(i=1;i<=2;i++)193.{194.for(ii=1;ii<=3;ii++)195.{196.h=3*(i-1)+ii;197.dh=3*(jm[e][i]-1)+ii;//给出整体坐标下单元节点位移198.wy[h]=p[dh];199.}200.}201.for(i=1;i<=6;i++)202.{203.f[i]=0.0;204.for(j=1;j<=6;j++)205.{206.for(k=1;k<=6;k++)//求由节点位移引起的单元节点力207.{208.f[i]=f[i]+kd[i][j]*t[j][k]*wy[k];209.}210.}211.}212.if(npf>0)213.{214.for(i=1;i<=npf;i++)//按非节点载荷数循环215.if(pf[i][3]==e)//找到荷载所在的单元216.{217.hz=i;218.gdnl(hz);//求固端反力219.for(j=1;j<=6;j++)//将固端反力累加220.{221.f[j]=f[j]+f0[j];222.}223.}224.}225.printf("%-3d(A)%9.5f%9.5f%9.5f\n",e,f[1],f[2],f[3]);//输出单元A （i）端力226.printf("(B)%9.5f%9.5f%9.5f\n",f[4],f[5],f[6]);//输出单元B（i）端力227.}228.return;229. }230. //**主程序结束**231.232. //************************************************************233. //<功能：将非节点载荷下的杆端力计算出来存入f0[]>234. //************************************************************235.236. void gdnl(int hz)237. {238.int ind,e;239.double g,c,l0,d;240.241.242.g=pf[hz][1];//载荷值243.c=pf[hz][2];//载荷位置244.e=(int)pf[hz][3];//作用单元245.ind=(int)pf[hz][4];//载荷类型246.l0=gc[e];//杆长247.d=l0-c;248.if(ind==1)249.{250.f0[1]=0.0;251.f0[2]=-(g*c*(2-2*c*c/(l0*l0)+(c*c*c)/(l0*l0*l0)))/2;//均布载荷的固端反力252.f0[3]=-(g*c*c)*(6-8*c/l0+3*c*c/(l0*l0))/12;253.f0[4]=0.0;254.f0[5]=-g*c-f0[2];255.f0[6]=(g*c*c*c)*(4-3*c/l0)/(12*l0);256.}257.else258.{259.if(ind==2)//横向集中力的固端反力260.{261.f0[1]=0.0;262.f0[2]=(-(g*d*d)*(l0+2*c))/(l0*l0*l0);263.f0[3]=-(g*c*d*d)/(l0*l0);264.f0[4]=0.0;265.f0[5]=(-g*c*c*(l0+2*d))/(l0*l0*l0);266.f0[6]=(g*c*c*d)/(l0*l0);267.}268.else269.{270.f0[1]=-(g*d/l0);//纵向集中力的固端反力271.f0[2]=0.0;272.f0[3]=0.0;273.f0[4]=-g*c/l0;274.f0[5]=0.0;275.f0[6]=0.0;276.}277.}278. }279.280. //************************************************************281. //<功能：构成坐标变换矩阵>282. //************************************************************283. void zb(int e)284. {285.double ceta,co,si;286.int i,j;287.ceta=(gj[e]*pi)/180;//角度变弧度288.co=cos(ceta);289.si=sin(ceta);290.t[1][1]=co;//计算T右上角元素291.t[1][2]=si;292.t[2][1]=-si;293.t[2][2]=co;294.t[3][3]=1.0;295.for(i=1;i<=3;i++)296.{297.for(j=1;j<=3;j++)//计算T的左下角元素298.{299.t[i+3][j+3]=t[i][j];300.}301.}302. }303.304.305.306. //*****************************************************307. //<计算局部坐标下单元刚度矩阵kd[][]>308. //*****************************************************309. void jdugd(int e)310. {311.double A0,l0,j0;312.int i;313.int j;314.315.316.A0=mj[e];//面积317.l0=gc[e];//杆长318.j0=gx[e];//惯性钜319.320.321.for(i=0;i<=6;i++)322.for(j=0;j<=6;j++)//kd清0323.kd[i][j]=0.0;324.325.kd[1][1]=e0*A0/l0;326.kd[2][2]=12*e0*j0/pow(l0,3);327.kd[3][2]=6*e0*j0/pow(l0,2);328.kd[3][3]=4*e0*j0/l0;329.kd[4][1]=-kd[1][1];330.kd[4][4]=kd[1][1];331.kd[5][2]=-kd[2][2];//计算kd左下角各元素332.kd[5][3]=-kd[3][2];333.kd[5][5]=kd[2][2];334.kd[6][2]=kd[3][2];335.kd[6][3]=2*e0*j0/l0;336.kd[6][5]=-kd[3][2];337.kd[6][6]=kd[3][3];338.339.for(i=1;i<=6;i++)340.for(j=1;j<=i;j++)//将kd左下角元素按对称原则送到右下角341.kd[j][i]=kd[i][j];342. }343.344.345. //**********************************************************346. //<计算整体坐标下单元刚度矩阵ke[][]>347. //**********************************************************348. void dugd(int e)349. {350.int i,k,j,m;351.jdugd(e);//计算局部单元刚度矩阵kd352.zb(e);//计算坐标变换矩阵T353.for(i=1;i<=6;i++)354.{355.for(j=1;j<=6;j++)356.{357.ke[i][j]=0.0;358.for(k=1;k<=6;k++)359.{360.for(m=1;m<=6;m++)361.{362.ke[i][j]=ke[i][j]+t[k][i]*kd[k][m]*t[m][j];//计算刚度坐标单元刚度矩阵ke 363.}364.}365.}366.}367. }368.369.370. //**程序结束**371.372.373.374. </math.h></stdio.h>。
有限元算例二维传热c++程序源代码4
//please turn to page 158 for more detailsprintf("计算基函数系数值...\n");for(id=0;id<nx*ny*2;id++){for(i=0;i<3;i++){if(i==0) j=1,k=2;else if(i==1) j=2,k=0;else if(i==2) j=0,k=1;pE[id].A=( (pE[id].nd[j].x-pE[id].nd[i].x)*(pE[id].nd[k].y-pE[id].nd[i].y)-(pE[id].nd[j].y-pE[id].nd[i].y)*(pE[id].nd[k].x-pE[id].nd[i].x) )/2.0;D=2.0*pE[id].A;pE[id].a[i]=( pE[id].nd[j].x*pE[id].nd[k].y- pE[id].nd[k].x*pE[id].nd[j].y )/D;pE[id].b[i]=( pE[id].nd[j].y-pE[id].nd[k].y )/D;pE[id].c[i]=( pE[id].nd[k].x-pE[id].nd[j].x )/D;}}printf("OK!\n");printf("计算单元有限元特征式系数矩阵...\n");int l,m;for(i=0;i<nx;i++) //计算单元有限元特征式系数矩阵for(j=0;j<ny;j++){for(l=0;l<3;l++) //for the first triangle in the rectanglefor(m=0;m<3;m++){pE[i*2+j*ny*2].Aij[l][m]=( pE[i*2+j*ny*2].b[l]*pE[i*2+j*ny*2].b[m] +pE[i*2+j*ny*2].c[l]*pE[i*2+j*ny*2].c[m] ) * pE[i*2+j*ny*2].A;}for(l=0;l<3;l++) //for the second triangle in the rectangle for(m=0;m<3;m++){pE[i*2+j*ny*2+1].Aij[l][m]=( pE[i*2+j*ny*2+1].b[l]*pE[i*2+j*ny*2+1].b[m] + pE[i*2+j*ny*2+1].c[l]*pE[i*2+j*ny*2+1].c[m] ) * pE[i*2+j*ny*2+1].A;}}printf("OK!\n");printf("计算积分值,填充到f函数向量数组...\n");static int js[2]={4,4}; //每一层积分区间均分为4个子区间int idx=0;for(i=0;i<nx;i++) //计算积分值,填充到f函数向量数组for(j=0;j<ny;j++){for(idx=0;idx<3;idx++) //for the first triangle in the rectangle。
有限元方法编程
有限元方法编程
【最新版】
目录
1.有限元方法概述
2.有限元方法的编程步骤
3.有限元方法的应用实例
4.总结
正文
一、有限元方法概述
有限元方法是一种数值分析方法，它通过将待求解的连续体划分为有限个小的、简单的子区域（单元），然后用这些单元的近似解描述整个连续体的行为。
这种方法主要用于求解偏微分方程，特别是在固体力学、流体力学、热传导等领域有着广泛的应用。
二、有限元方法的编程步骤
1.几何建模：首先需要对问题进行几何建模，即将实际问题转化为计算机可以处理的数学模型。
这包括对物体的边界、形状等进行描述。
2.网格划分：将整个模型划分为有限个小的单元，这些单元可以是四面体、六面体等，根据问题的实际情况和求解的需要来选择。
3.选择适当的有限元公式：根据问题的性质和求解的目标，选择合适的有限元公式来描述单元内的物理量，如应力、应变等。
4.组装方程：将所有单元的公式组合起来，得到整个模型的方程。
5.求解方程：通过数值方法（如迭代法）求解得到的方程组，得到模型的解。
6.后处理：对求解结果进行分析和处理，如绘制应力分布图、应变分
布图等。
三、有限元方法的应用实例
有限元方法在许多工程领域都有广泛的应用，如飞机设计、桥梁设计、汽车设计等。
例如，在飞机设计中，可以通过有限元方法求解机翼的应力分布，从而优化机翼的设计，提高飞行性能。
四、总结
有限元方法是一种强大的数值分析工具，它可以用于求解各种复杂的工程问题。
通过几何建模、网格划分、选择适当的有限元公式、组装方程、求解方程和后处理等步骤，可以得到问题的解。
三结点三角形有限单元程序设计C 语言
fclose(fp) ;
int NN2=NN*2; int NBW=NN2;
//半带宽 NBW
float GK[50][50];
float DU[50];
int *count=new int [NN];
float (*strain)[3]=new float [NN][3];
float (*stress)[3]=new float [NN][3];
int ngn[3]; for(i=0;i<3;i++)ngn[i]=NGN[NOE-1][i];
float cn[3][2]; for(i=0;i<3;i++) for(j=0;j<2;j++) cn[i][j]=CN[ngn[i]-1][j];
float a[3]={0,0,0}; float b[3]={0,0,0}; float c[3]={0,0,0};
ｖｏｉｄ ｅｌｅｓｔｒａｉｎ＿ｓｔｒｅｓｓ（
ｉｎｔＮＯＥ，ｉｎｔＮＧＮ［］［３］，ｆｌｏａｔＣＮ［］［２］，
ｆｌｏａｔＤＵ［５０］，ｆｌｏａｔＥ０，ｆｌｏａｔＵ０，ｆｌｏａｔｔ，
ｉｎｔＮＥ，ｉｎｔＮＮ，ｆｌｏａｔｅｌｅｓｔｒａｉｎ［３］，ｆｌｏａｔｅｌｅｓｔｒｅｓｓ［３］）；
＃ｅｎｄｉｆ
Ⅴ．节点应力－应变头文件 ＮｏｄｅＳｔｒａｉｎＳｔｒｅｓｓ．ｈ
fscanf(fp,"%s",&newlin); for(i=0;i<NN;i++)
for(j=0;j<2;j++) fscanf(fp, "%f", &CN[s",&newlin); int NRE=0; fscanf(fp, "%d", &NRE) ; int *restrict=new int [NRE]; for(i=0;i<NRE;i++)
Matlab 有限元法计算分析程序编写
MATLAB的使用方法
1) 最简单的计算器使用法 求[12＋2×(7－4)]÷32的算术运算结果 （1）用键盘在MATLAB指令窗中输入一下内容 (12+2*(7-4))/3^2 （2）在上述表达式输入完成后，按【Enter】键，该指令被执行 （3）在指令执行后，MATLAB指令窗中将显示一下内容 ans = 2 [说明] 加 + 减 乘 * 除 / 或 \ (这两个符号对于数组有不同的含义) 幂 ^ “ans”是answer的缩写，其含义是运算答案，它是MATLAB的一个默 认变量
material_number = fscanf( fid_in, '%d', 1 ) ; % read material number material = zeros( material_number, 3 ) ; for i=1:1:material_number nm = fscanf( fid_in, '%d', 1 ) ; material( i, : ) = fscanf( fid_in, '%f', [1,3] ) ; % read materials definition end bc_number = fscanf( fid_in, '%d', 1 ) ; % read boundary conditions number bc = zeros( bc_number, 3 ) ; for i=1:1:bc_number % read boundary condition definition bc( i, 1 ) = fscanf( fid_in, '%d', 1 ) ; bc( i, 2 ) = fscanf( fid_in, '%d', 1 ) ; bc( i, 3 ) = fscanf( fid_in, '%f', 1 ) ; end
有限元分析中《结构力学》矩阵位移法C语言程序(附例题)
程序：#include "stdafx.h"#include "stdio.h"#include "math.h"#include "stdlib.h"void main(){int loc[3][2]={0},ifix[6]={0};float area[3]={0.0},fint[3]={0.0},cx[4]={0.0},cy[4]={0.0},f[12]={0.0},fr[12]= {0.0},fe[3][6]={0.0};int nn,ne,nd,nfix;float ea;int i,j,k;FILE *shuru,*shuchu;shuru=fopen("shuru.dat","r");shuchu=fopen("shuchu.dat","w");fscanf(shuru,"%d%d%d%d%f",&nn,&ne,&nd,&nfix,&ea);fprintf(shuchu,"nn ne nd nfix e\n%d %d %d %d %f\n",nn,ne,nd,nfix,ea);i=0;while(i<=ne-1){fscanf(shuru,"%d%d%f%f",&loc[i][0],&loc[i][1],&area[i],&fint[i]);i++;}fprintf(shuchu,"element node1 node2 area fint\n");i=0;while(i<=ne-1){fprintf(shuchu,"%d %d %d %f %f\n",i+1,loc[i][0],loc[i][1],area[i],fint[i]);i++;}j=0;while(j<=nn-1){fscanf(shuru,"%f%f",&cx[j],&cy[j]);j++;}fprintf(shuchu,"node x-coord y-coord\n");j=0;while(j<=nn-1){fprintf(shuchu,"%d %f %f\n",j+1,cx[j],cy[j]);j++;}k=0;while(k<=nfix-1){fscanf(shuru,"%d",&ifix[k]);k++;}fprintf(shuchu,"ifix=");k=0;while(k<=nfix-1){fprintf(shuchu,"%d ",ifix[k]);k++;}fprintf(shuchu,"\n");void cst(int (*loc)[2],int *ifix,float *area,float *fint,float *cx,float *cy,float*f,float *fr,float (*fe)[6],FILE *shuru,FILE *shuchu,float ea);cst(loc,ifix,area,fint,cx,cy,f,fr,fe,shuru,shuchu,ea);fprintf(shuchu,"node x-disp y-disp thita\n");i=0;while(i<=3){fprintf(shuchu,"%d %f %f %f\n",i+1,f[3*i],f[3*i+1],f[3*i+2]);i++;}fprintf(shuchu,"reaction nodal forces from the equations\n");fprintf(shuchu,"node x-load y-load moment\n");i=0;while(i<=3){fprintf(shuchu,"%d %f %f %f\n",i+1,fr[3*i],fr[3*i+1],fr[3*i+2]);i++;}fprintf(shuchu,"element axi-f shear-q moment-m\n");i=0;while(i<=ne-1){fprintf(shuchu,"%d %f %f %f %f %f %f\n",i+1,fe[i][0],fe[i][1],fe[i][2],fe[i][3],fe [i][4],fe[i][5]);i++;}fclose(shuru);fclose(shuchu);}void cst(int (*loc)[2],int *ifix,float *area,float *fint,float *cx,float *cy,float*f,float *fr,float (*fe)[6],FILE *shuru,FILE *shuchu,float ea){int np,nvd;float p1[3][6]={0.0},p2[3][6]={0.0},gk[12][12]={0.0},gk1[12][12]={0.0},al[3]= {0.0},tt[3][6][6]={0.0},bkl[3][6][6]={0.0};float t[6][6]={0.0},css[3]={0.0},snn[3]={0.0},ek[6][6]={0.0},ekl[6][6]={0.0},ekk[3] [6][6]={0.0},xx[6]={0.0},ba[6][6]={0.0};int nn=4,ne=3,nd=12,nfix=6;int ii,jj,i,j,k,l,inode,nodei,idofn,nrows,nrowe,jnode,nodej,jdofn,ncols,ncole;int i1,i2,ie,ix;float x12,y12,q,eal,eil1,eil2,eil3;i=0;{for(;i<=ne-1;i++){i1=loc[i][0];i2=loc[i][1];x12=cx[i2-1]-cx[i1-1];y12=cy[i2-1]-cy[i1-1];al[i]=sqrt(pow(x12,2)+pow(y12,2));css[i]=x12/al[i];snn[i]=y12/al[i];}}fscanf(shuru,"%d%d",&np,&nvd);if(np!=0)i=0;for(;i<=np-1;i++){fscanf(shuru,"%d%f%f%f",&i,&f[3*i],&f[3*i+1],&f[3*i+2]); }if(nvd!=0)i=0;for(;i<=nvd-1;i++){fscanf(shuru,"%d%f",&ie,&q);i1=loc[ie-1][0];i2=loc[ie-1][1];p1[ie-1][1]=q*al[ie-1]/2;p1[ie-1][2]=q*al[ie-1]*al[ie-1]/12;p1[ie-1][4]=q*al[ie-1]/2;p1[ie-1][5]=-q*al[ie-1]*al[ie-1]/12;}i=0;for(;i<=nd-1;i++){j=0;for(;j<=nd-1;j++){gk[i][j]=0.0;}}for(i=0;i<=ne-1;i++){j=0;for(;j<=5;j++){k=0;for(;k<=5;k++){ekl[j][k]=0.0;ek[j][k]=0.0;t[j][k]=0.0;}}eal=ea*area[i]/al[i];eil1=ea*fint[i]/al[i];eil2=ea*fint[i]/(al[i]*al[i]);eil3=ea*fint[i]/(al[i]*al[i]*al[i]);ekl[0][0]=eal;ekl[1][1]=12.0*eil3;ekl[2][2]=4.0*eil1;ekl[3][3]=eal;ekl[4][4]=12.0*eil3;ekl[5][5]=4.0*eil1;ekl[2][1]=6.0*eil2;ekl[3][0]=-eal;ekl[4][1]=-12.0*eil3;ekl[4][2]=-6.0*eil2;ekl[5][1]=6.0*eil2;ekl[5][2]=2.0*eil1;ekl[5][4]=-6.0*eil2;ii=0;for(;ii<=4;ii++){jj=ii+1;for(;jj<=5;jj++){ekl[ii][jj]=ekl[jj][ii];}}k=0;for(;k<=5;k++){l=0;for(;l<=5;l++){{ekk[i][k][l]=ekl[k][l];fprintf(shuchu,"%d %d %d %f %f\n",i+1,k+1,l+1,ekl[k][l],ekk[i][k][l]);} }}t[0][0]=css[i];t[0][1]=-snn[i];t[1][0]=snn[i];t[1][1]=css[i];t[2][2]=1.0;t[3][3]=css[i];t[3][4]=-snn[i];t[4][3]=snn[i];t[4][4]=css[i];t[5][5]=1.0;j=0;for(;j<=5;j++){k=0;for(;k<=5;k++){tt[i][j][k]=t[j][k];p2[i][j]=p2[i][j]+t[j][k]*p1[i][k];}}ii=0;for(;ii<=5;ii++){j=0;for(;j<=5;j++){ba[ii][j]=0.0;k=0;for(;k<=5;k++){ba[ii][j]=ba[ii][j]+tt[i][ii][k]*ekl[k][j];}}}ek[ii][j]=0.0;ii=0;for(;ii<=5;ii++){j=0;for(;j<=5;j++){k=0;for(;k<=5;k++){ek[ii][j]=ek[ii][j]+ba[ii][k]*tt[i][j][k];}}}j=0;for(;j<=5;j++){ii=0;for(;ii<=5;ii++){fprintf(shuchu,"i,ii,j,ek,tt=%d %d %d %f %f\n",i+1,ii+1,j+1,ek[ii][j],tt[i][ii][j]); }}inode=0;while(inode<=1){nodei=loc[i][inode];idofn=0;while(idofn<=2){nrows=(nodei-1)*3+idofn;nrowe=inode*3+idofn;f[nrows]=f[nrows]+p2[i][nrowe];jnode=0;while(jnode<=1){nodej=loc[i][jnode];jdofn=0;while(jdofn<=2){ncols=(nodej-1)*3+jdofn;ncole=jnode*3+jdofn;gk[nrows][ncols]=gk[nrows][ncols]+ek[nrowe][ncole];jdofn++;}jnode++;}idofn++;}inode++;}}i=0;for(;i<=nd-1;i++){j=0;for(;j<=nd-1;j++){gk1[i][j]=gk[i][j];}}fprintf(shuchu,"nodal forces from applied loads\node x-load y-load moment\n"); i=0;for(;i<=nn-1;i++){fprintf(shuchu,"%d %f %f %f\n",i+1,f[3*i],f[3*i+1],f[3*i+2]);}i=0;for(;i<=nd-1;i++){j=0;for(;j<=nd-1;j++){fprintf(shuchu,"i,j,gk1 %d %d %f\n",i+1,j+1,gk[i][j]);}}i=0;for(;i<=nfix-1;i++){ix=ifix[i];gk[ix-1][ix-1]=gk[ix-1][ix-1]*1.0e20;}void gauss(float (*a)[12],float *b,int n);gauss(gk,f,nd);i=0;for(;i<=nd-1;i++){fr[i]=0.0;j=0;for(;j<=nd-1;j++){fr[i]=fr[i]+gk1[i][j]*f[j];fr[i]=fr[i]-f[i];}}for(i=0;i<=ne-1;i++){i1=loc[i][0];i2=loc[i][1];xx[0]=f[3*i1-3];xx[1]=f[3*i1-2];xx[2]=f[3*i1-1];xx[3]=f[3*i2-3];xx[4]=f[3*i2-2];xx[5]=f[3*i2-1];j=0;for(;j<=5;j++){k=0;for(;k<=5;k++){bkl[i][j][k]=0.0;l=0;for(;l<=5;l++){bkl[i][j][k]=bkl[i][j][k]+ekk[i][j][l]*tt[i][k][l];}}}j=0;for(;j<=5;j++){fe[i][j]=0.0;k=0;for(;k<=5;k++){fe[i][j]=fe[i][j]+bkl[i][j][k]*xx[k];}}j=0;for(;j<=5;j++){fe[i][j]=fe[i][j]-p1[i][j];}}}void gauss(float (*a)[12],float *b,int n){int i,i1,j,m;i=0;for(;i<=n-1;i++){i1=i+1;for(j=i1;j<=n-1;j++)a[i][j]=a[i][j]/a[i][i];b[i]=b[i]/a[i][i];a[i][i]=1.0;for(j=i1;j<=n-1;j++){for(m=i1;m<=n-1;m++)a[j][m]=a[j][m]-a[j][i]*a[i][m];b[j]=b[j]-a[j][i]*b[i];}}i=n-2;for(;i>=0;i--){j=i+1;for(;j<=n-1;j++){b[i]=b[i]-a[i][j]*b[j];}}}一、如图所示平面刚架的内力，各杆面积A=76.3cm2,惯性矩I=15760cm4，弹性模量E=2×105MPa程序：#include "stdafx.h"#include "stdio.h"#include "math.h"#include "stdlib.h"void main(){int loc[3][2]={0},ifix[6]={0};float area[3]={0.0},fint[3]={0.0},cx[4]={0.0},cy[4]={0.0},f[12]={0.0},fr[12]={0.0},fe[3][6]={0.0};int nn,ne,nd,nfix;float ea;int i,j,k;FILE *shuru,*shuchu;shuru=fopen("shuru.dat","r");shuchu=fopen("shuchu.dat","w");fscanf(shuru,"%d%d%d%d%f",&nn,&ne,&nd,&nfix,&ea);fprintf(shuchu,"nn ne nd nfix e\n%d %d %d %d %f\n",nn,ne,nd,nfix,ea);i=0;while(i<=ne-1){fscanf(shuru,"%d%d%f%f",&loc[i][0],&loc[i][1],&area[i],&fint[i]);fprintf(shuchu,"element node1 node2 area fint\n");i=0;while(i<=ne-1){fprintf(shuchu,"%d %d %d %f %f\n",i+1,loc[i][0],loc[i][1],area[i],fint[i]);i++;}j=0;while(j<=nn-1){fscanf(shuru,"%f%f",&cx[j],&cy[j]);j++;}fprintf(shuchu,"node x-coord y-coord\n");j=0;while(j<=nn-1){fprintf(shuchu,"%d %f %f\n",j+1,cx[j],cy[j]);j++;}k=0;while(k<=nfix-1){fscanf(shuru,"%d",&ifix[k]);k++;}fprintf(shuchu,"ifix=");k=0;while(k<=nfix-1){fprintf(shuchu,"%d ",ifix[k]);k++;}fprintf(shuchu,"\n");void cst(int (*loc)[2],int *ifix,float *area,float *fint,float *cx,float *cy,float*f,float *fr,float (*fe)[6],FILE *shuru,FILE *shuchu,float ea);cst(loc,ifix,area,fint,cx,cy,f,fr,fe,shuru,shuchu,ea);fprintf(shuchu,"node x-disp y-disp thita\n");i=0;while(i<=3){fprintf(shuchu,"%d %f %f %f\n",i+1,f[3*i],f[3*i+1],f[3*i+2]);i++;}fprintf(shuchu,"reaction nodal forces from the equations\n");fprintf(shuchu,"node x-load y-load moment\n");i=0;while(i<=3){fprintf(shuchu,"%d %f %f %f\n",i+1,fr[3*i],fr[3*i+1],fr[3*i+2]);i++;}fprintf(shuchu,"element axi-f shear-q moment-m\n");i=0;while(i<=ne-1){fprintf(shuchu,"%d %f %f %f %f %f %f\n",i+1,fe[i][0],fe[i][1],fe[i][2],fe[i][3],fe [i][4],fe[i][5]);fclose(shuru);fclose(shuchu);}void cst(int (*loc)[2],int *ifix,float *area,float *fint,float *cx,float *cy,float*f,float *fr,float (*fe)[6],FILE *shuru,FILE *shuchu,float ea){int np,nvd;float p1[3][6]={0.0},p2[3][6]={0.0},gk[12][12]={0.0},gk1[12][12]={0.0},al[3]= {0.0},tt[3][6][6]={0.0},bkl[3][6][6]={0.0};float t[6][6]={0.0},css[3]={0.0},snn[3]={0.0},ek[6][6]={0.0},ekl[6][6]={0.0},ekk[3] [6][6]={0.0},xx[6]={0.0},ba[6][6]={0.0};int nn=4,ne=3,nd=12,nfix=6;int ii,jj,i,j,k,l,inode,nodei,idofn,nrows,nrowe,jnode,nodej,jdofn,ncols,ncole;int i1,i2,ie,ix;float x12,y12,q,eal,eil1,eil2,eil3;i=0;{for(;i<=ne-1;i++){i1=loc[i][0];i2=loc[i][1];x12=cx[i2-1]-cx[i1-1];y12=cy[i2-1]-cy[i1-1];al[i]=sqrt(pow(x12,2)+pow(y12,2));css[i]=x12/al[i];snn[i]=y12/al[i];}}fscanf(shuru,"%d%d",&np,&nvd);if(np!=0)i=0;for(;i<=np-1;i++){fscanf(shuru,"%d%f%f%f",&i,&f[3*i],&f[3*i+1],&f[3*i+2]);}if(nvd!=0)i=0;for(;i<=nvd-1;i++){fscanf(shuru,"%d%f",&ie,&q);i1=loc[ie-1][0];i2=loc[ie-1][1];p1[ie-1][1]=q*al[ie-1]/2;p1[ie-1][2]=q*al[ie-1]*al[ie-1]/12;p1[ie-1][4]=q*al[ie-1]/2;p1[ie-1][5]=-q*al[ie-1]*al[ie-1]/12;}i=0;for(;i<=nd-1;i++){j=0;for(;j<=nd-1;j++){gk[i][j]=0.0;}}for(i=0;i<=ne-1;i++){j=0;for(;j<=5;j++){k=0;for(;k<=5;k++){ekl[j][k]=0.0;ek[j][k]=0.0;t[j][k]=0.0;}}eal=ea*area[i]/al[i];eil1=ea*fint[i]/al[i];eil2=ea*fint[i]/(al[i]*al[i]);eil3=ea*fint[i]/(al[i]*al[i]*al[i]); ekl[0][0]=eal;ekl[1][1]=12.0*eil3;ekl[2][2]=4.0*eil1;ekl[3][3]=eal;ekl[4][4]=12.0*eil3;ekl[5][5]=4.0*eil1;ekl[2][1]=6.0*eil2;ekl[3][0]=-eal;ekl[4][1]=-12.0*eil3;ekl[4][2]=-6.0*eil2;ekl[5][1]=6.0*eil2;ekl[5][2]=2.0*eil1;ekl[5][4]=-6.0*eil2;ii=0;for(;ii<=4;ii++){jj=ii+1;for(;jj<=5;jj++){ekl[ii][jj]=ekl[jj][ii];}}k=0;for(;k<=5;k++){l=0;for(;l<=5;l++){{ekk[i][k][l]=ekl[k][l];fprintf(shuchu,"%d %d %d %f %f\n",i+1,k+1,l+1,ekl[k][l],ekk[i][k][l]);} }}t[0][0]=css[i];t[0][1]=-snn[i];t[1][0]=snn[i];t[1][1]=css[i];t[2][2]=1.0;t[3][3]=css[i];t[3][4]=-snn[i];t[4][3]=snn[i];t[4][4]=css[i];t[5][5]=1.0;j=0;for(;j<=5;j++){k=0;for(;k<=5;k++){tt[i][j][k]=t[j][k];p2[i][j]=p2[i][j]+t[j][k]*p1[i][k];}}ii=0;for(;ii<=5;ii++){j=0;for(;j<=5;j++){ba[ii][j]=0.0;k=0;for(;k<=5;k++){ba[ii][j]=ba[ii][j]+tt[i][ii][k]*ekl[k][j];}}}ek[ii][j]=0.0;ii=0;for(;ii<=5;ii++){j=0;for(;j<=5;j++){k=0;for(;k<=5;k++){ek[ii][j]=ek[ii][j]+ba[ii][k]*tt[i][j][k];}}}j=0;for(;j<=5;j++){ii=0;for(;ii<=5;ii++){fprintf(shuchu,"i,ii,j,ek,tt=%d %d %d %f %f\n",i+1,ii+1,j+1,ek[ii][j],tt[i][ii][j]); }}inode=0;while(inode<=1){nodei=loc[i][inode];idofn=0;while(idofn<=2){nrows=(nodei-1)*3+idofn;nrowe=inode*3+idofn;f[nrows]=f[nrows]+p2[i][nrowe];jnode=0;while(jnode<=1){nodej=loc[i][jnode];jdofn=0;while(jdofn<=2){ncols=(nodej-1)*3+jdofn;ncole=jnode*3+jdofn;gk[nrows][ncols]=gk[nrows][ncols]+ek[nrowe][ncole];jdofn++;}jnode++;}idofn++;}inode++;}}i=0;for(;i<=nd-1;i++){j=0;for(;j<=nd-1;j++){gk1[i][j]=gk[i][j];}}fprintf(shuchu,"nodal forces from applied loads\node x-load y-load moment\n"); i=0;for(;i<=nn-1;i++){fprintf(shuchu,"%d %f %f %f\n",i+1,f[3*i],f[3*i+1],f[3*i+2]);}i=0;for(;i<=nd-1;i++){j=0;for(;j<=nd-1;j++){fprintf(shuchu,"i,j,gk1 %d %d %f\n",i+1,j+1,gk[i][j]); }}i=0;for(;i<=nfix-1;i++){ix=ifix[i];gk[ix-1][ix-1]=gk[ix-1][ix-1]*1.0e20;}void gauss(float (*a)[12],float *b,int n);gauss(gk,f,nd);i=0;for(;i<=nd-1;i++){fr[i]=0.0;j=0;for(;j<=nd-1;j++){fr[i]=fr[i]+gk1[i][j]*f[j];fr[i]=fr[i]-f[i];}}for(i=0;i<=ne-1;i++){i1=loc[i][0];i2=loc[i][1];xx[0]=f[3*i1-3];xx[1]=f[3*i1-2];xx[2]=f[3*i1-1];xx[3]=f[3*i2-3];xx[4]=f[3*i2-2];xx[5]=f[3*i2-1];j=0;for(;j<=5;j++){k=0;for(;k<=5;k++){bkl[i][j][k]=0.0;l=0;for(;l<=5;l++){bkl[i][j][k]=bkl[i][j][k]+ekk[i][j][l]*tt[i][k][l];}}}j=0;for(;j<=5;j++){fe[i][j]=0.0;k=0;for(;k<=5;k++){fe[i][j]=fe[i][j]+bkl[i][j][k]*xx[k]; }}j=0;for(;j<=5;j++){fe[i][j]=fe[i][j]-p1[i][j];}}}void gauss(float (*a)[12],float *b,int n) {int i,i1,j,m;i=0;for(;i<=n-1;i++){i1=i+1;for(j=i1;j<=n-1;j++)a[i][j]=a[i][j]/a[i][i];b[i]=b[i]/a[i][i];a[i][i]=1.0;for(j=i1;j<=n-1;j++){for(m=i1;m<=n-1;m++)a[j][m]=a[j][m]-a[j][i]*a[i][m];b[j]=b[j]-a[j][i]*b[i];}}i=n-2;for(;i>=0;i--){j=i+1;for(;j<=n-1;j++){b[i]=b[i]-a[i][j]*b[j];}}}程序输出：nn ne nd nfix e4 3 12 6 200000000.000000element node1 node2 area fint1 12 0.007630 0.0001582 3 1 0.007630 0.0001583 24 0.007630 0.000158node x-coord y-coord1 0.000000 5.0000002 6.400000 5.0000003 0.000000 0.0000004 9.600000 0.000000ifix=7 8 9 10 11 121 1 1 238437.500000 238437.5000001 12 0.000000 0.0000001 1 3 0.000000 0.0000001 1 4 -238437.500000 -238437.5000001 1 5 0.000000 0.0000001 1 6 0.000000 0.0000001 2 1 0.000000 0.0000001 2 2 1442.871094 1442.8710941 2 3 4617.187500 4617.1875001 2 4 0.000000 0.0000001 2 5 -1442.871094 -1442.8710941 2 6 4617.187500 4617.1875001 3 1 0.000000 0.0000001 32 4617.187500 4617.1875001 3 3 19700.000000 19700.0000001 3 4 0.000000 0.0000001 3 5 -4617.187500 -4617.1875001 3 6 9850.000000 9850.0000001 4 1 -238437.500000 -238437.5000001 42 0.000000 0.0000001 4 3 0.000000 0.0000001 4 4 238437.500000 238437.5000001 4 5 0.000000 0.0000001 4 6 0.000000 0.0000001 5 1 0.000000 0.0000001 52 -1442.871094 -1442.8710941 5 3 -4617.187500 -4617.1875001 5 4 0.000000 0.0000001 5 5 1442.871094 1442.8710941 5 6 -4617.187500 -4617.1875001 6 1 0.000000 0.0000001 62 4617.187500 4617.1875001 6 3 9850.000000 9850.0000001 6 4 0.000000 0.0000001 6 5 -4617.187500 -4617.1875001 6 6 19700.000000 19700.000000i,ii,j,ek,tt=1 1 1 238437.500000 1.000000 i,ii,j,ek,tt=1 2 1 0.000000 0.000000i,ii,j,ek,tt=1 3 1 0.000000 0.000000i,ii,j,ek,tt=1 4 1 -238437.500000 0.000000 i,ii,j,ek,tt=1 5 1 0.000000 0.000000i,ii,j,ek,tt=1 6 1 0.000000 0.000000i,ii,j,ek,tt=1 1 2 0.000000 0.000000i,ii,j,ek,tt=1 2 2 1442.871094 1.000000i,ii,j,ek,tt=1 3 2 4617.187500 0.000000 i,ii,j,ek,tt=1 4 2 0.000000 0.000000i,ii,j,ek,tt=1 5 2 -1442.871094 0.000000 i,ii,j,ek,tt=1 6 2 4617.187500 0.000000 i,ii,j,ek,tt=1 1 3 0.000000 0.000000i,ii,j,ek,tt=1 2 3 4617.187500 0.000000 i,ii,j,ek,tt=1 3 3 19700.000000 1.000000 i,ii,j,ek,tt=1 4 3 0.000000 0.000000i,ii,j,ek,tt=1 5 3 -4617.187500 0.000000 i,ii,j,ek,tt=1 6 3 9850.000000 0.000000 i,ii,j,ek,tt=1 1 4 -238437.500000 0.000000 i,ii,j,ek,tt=1 2 4 0.000000 0.000000i,ii,j,ek,tt=1 3 4 0.000000 0.000000i,ii,j,ek,tt=1 4 4 238437.500000 1.000000 i,ii,j,ek,tt=1 5 4 0.000000 0.000000i,ii,j,ek,tt=1 6 4 0.000000 0.000000i,ii,j,ek,tt=1 1 5 0.000000 0.000000i,ii,j,ek,tt=1 2 5 -1442.871094 0.000000 i,ii,j,ek,tt=1 3 5 -4617.187500 0.000000 i,ii,j,ek,tt=1 4 5 0.000000 0.000000i,ii,j,ek,tt=1 5 5 1442.871094 1.000000 i,ii,j,ek,tt=1 6 5 -4617.187500 0.000000 i,ii,j,ek,tt=1 1 6 0.000000 0.000000i,ii,j,ek,tt=1 2 6 4617.187500 0.000000 i,ii,j,ek,tt=1 3 6 9850.000000 0.000000 i,ii,j,ek,tt=1 4 6 0.000000 0.000000i,ii,j,ek,tt=1 5 6 -4617.187500 0.000000 i,ii,j,ek,tt=1 6 6 19700.000000 1.000000 2 1 1 305200.000000 305200.0000002 1 2 0.000000 0.0000002 13 0.000000 0.0000002 1 4 -305200.000000 -305200.0000002 1 5 0.000000 0.0000002 1 6 0.000000 0.0000002 2 1 0.000000 0.0000002 2 2 3025.919922 3025.9199222 23 7564.800293 7564.8002932 2 4 0.000000 0.0000002 2 5 -3025.919922 -3025.9199222 2 6 7564.800293 7564.8002932 3 1 0.000000 0.0000002 3 2 7564.800293 7564.8002932 3 3 25216.001953 25216.0019532 3 4 0.000000 0.0000002 3 5 -7564.800293 -7564.8002932 3 6 12608.000977 12608.0009772 4 1 -305200.000000 -305200.0000002 4 2 0.000000 0.0000002 43 0.000000 0.0000002 4 4 305200.000000 305200.0000002 4 5 0.000000 0.0000002 4 6 0.000000 0.0000002 5 1 0.000000 0.0000002 5 2 -3025.919922 -3025.9199222 53 -7564.800293 -7564.8002932 5 4 0.000000 0.0000002 5 5 3025.919922 3025.9199222 5 6 -7564.800293 -7564.8002932 6 1 0.000000 0.0000002 6 2 7564.800293 7564.8002932 63 12608.000977 12608.0009772 6 4 0.000000 0.0000002 6 5 -7564.800293 -7564.8002932 6 6 25216.001953 25216.001953i,ii,j,ek,tt=2 1 1 3025.919922 0.000000 i,ii,j,ek,tt=2 2 1 0.000000 1.000000i,ii,j,ek,tt=2 3 1 -7564.800293 0.000000 i,ii,j,ek,tt=2 4 1 -3025.919922 0.000000 i,ii,j,ek,tt=2 5 1 0.000000 0.000000i,ii,j,ek,tt=2 6 1 -7564.800293 0.000000 i,ii,j,ek,tt=2 1 2 0.000000 -1.000000i,ii,j,ek,tt=2 2 2 305200.000000 0.000000 i,ii,j,ek,tt=2 3 2 0.000000 0.000000i,ii,j,ek,tt=2 4 2 0.000000 0.000000i,ii,j,ek,tt=2 5 2 -305200.000000 0.000000 i,ii,j,ek,tt=2 6 2 0.000000 0.000000i,ii,j,ek,tt=2 1 3 -7564.800293 0.000000 i,ii,j,ek,tt=2 2 3 0.000000 0.000000i,ii,j,ek,tt=2 3 3 25216.001953 1.000000 i,ii,j,ek,tt=2 4 3 7564.800293 0.000000 i,ii,j,ek,tt=2 5 3 0.000000 0.000000i,ii,j,ek,tt=2 6 3 12608.000977 0.000000 i,ii,j,ek,tt=2 1 4 -3025.919922 0.000000 i,ii,j,ek,tt=2 2 4 0.000000 0.000000i,ii,j,ek,tt=2 3 4 7564.800293 0.000000 i,ii,j,ek,tt=2 4 4 3025.919922 0.000000 i,ii,j,ek,tt=2 5 4 0.000000 1.000000i,ii,j,ek,tt=2 6 4 7564.800293 0.000000i,ii,j,ek,tt=2 1 5 0.000000 0.000000i,ii,j,ek,tt=2 2 5 -305200.000000 0.000000 i,ii,j,ek,tt=2 3 5 0.000000 0.000000i,ii,j,ek,tt=2 4 5 0.000000 -1.000000i,ii,j,ek,tt=2 5 5 305200.000000 0.000000 i,ii,j,ek,tt=2 6 5 0.000000 0.000000i,ii,j,ek,tt=2 1 6 -7564.800293 0.000000 i,ii,j,ek,tt=2 2 6 0.000000 0.000000i,ii,j,ek,tt=2 3 6 12608.000977 0.000000 i,ii,j,ek,tt=2 4 6 7564.800293 0.000000 i,ii,j,ek,tt=2 5 6 0.000000 0.000000i,ii,j,ek,tt=2 6 6 25216.001953 1.000000 3 1 1 257061.218750 257061.2187503 1 2 0.000000 0.0000003 1 3 0.000000 0.0000003 14 -257061.218750 -257061.2187503 1 5 0.000000 0.0000003 1 6 0.000000 0.0000003 2 1 0.000000 0.0000003 2 2 1808.063232 1808.0632323 2 3 5366.628906 5366.6289063 24 0.000000 0.0000003 2 5 -1808.063232 -1808.0632323 2 6 5366.628906 5366.6289063 3 1 0.000000 0.0000003 3 2 5366.628906 5366.6289063 3 3 21238.716797 21238.7167973 34 0.000000 0.0000003 3 5 -5366.628906 -5366.6289063 3 6 10619.358398 10619.3583983 4 1 -257061.218750 -257061.2187503 4 2 0.000000 0.0000003 4 3 0.000000 0.0000003 4 4 257061.218750 257061.2187503 4 5 0.000000 0.0000003 4 6 0.000000 0.0000003 5 1 0.000000 0.0000003 5 2 -1808.063232 -1808.0632323 5 3 -5366.628906 -5366.6289063 54 0.000000 0.0000003 5 5 1808.063232 1808.0632323 5 6 -5366.628906 -5366.6289063 6 1 0.000000 0.0000003 6 2 5366.628906 5366.6289063 6 3 10619.358398 10619.3583983 64 0.000000 0.0000003 6 5 -5366.628906 -5366.6289063 6 6 21238.716797 21238.716797i,ii,j,ek,tt=3 1 1 75979.257813 0.539054 i,ii,j,ek,tt=3 2 1 -115892.476563 -0.842271 i,ii,j,ek,tt=3 3 1 4520.158203 0.000000i,ii,j,ek,tt=3 4 1 -75979.257813 0.000000 i,ii,j,ek,tt=3 5 1 115892.476563 0.000000 i,ii,j,ek,tt=3 6 1 4520.158203 0.000000i,ii,j,ek,tt=3 1 2 -115892.476563 0.842271 i,ii,j,ek,tt=3 2 2 182890.046875 0.539054 i,ii,j,ek,tt=3 3 2 2892.901367 0.000000i,ii,j,ek,tt=3 4 2 115892.476563 0.000000 i,ii,j,ek,tt=3 5 2 -182890.046875 0.000000 i,ii,j,ek,tt=3 6 2 2892.901367 0.000000i,ii,j,ek,tt=3 1 3 4520.158203 0.000000i,ii,j,ek,tt=3 2 3 2892.901367 0.000000i,ii,j,ek,tt=3 3 3 21238.716797 1.000000 i,ii,j,ek,tt=3 4 3 -4520.158203 0.000000 i,ii,j,ek,tt=3 5 3 -2892.901367 0.000000 i,ii,j,ek,tt=3 6 3 10619.358398 0.000000 i,ii,j,ek,tt=3 1 4 -75979.257813 0.000000 i,ii,j,ek,tt=3 2 4 115892.476563 0.000000 i,ii,j,ek,tt=3 3 4 -4520.158203 0.000000 i,ii,j,ek,tt=3 4 4 75979.257813 0.539054 i,ii,j,ek,tt=3 5 4 -115892.476563 -0.842271 i,ii,j,ek,tt=3 6 4 -4520.158203 0.000000 i,ii,j,ek,tt=3 1 5 115892.476563 0.000000 i,ii,j,ek,tt=3 2 5 -182890.046875 0.000000 i,ii,j,ek,tt=3 3 5 -2892.901367 0.000000 i,ii,j,ek,tt=3 4 5 -115892.476563 0.842271 i,ii,j,ek,tt=3 5 5 182890.046875 0.539054 i,ii,j,ek,tt=3 6 5 -2892.901367 0.000000 i,ii,j,ek,tt=3 1 6 4520.158203 0.000000i,ii,j,ek,tt=3 2 6 2892.901367 0.000000i,ii,j,ek,tt=3 3 6 10619.358398 0.000000 i,ii,j,ek,tt=3 4 6 -4520.158203 0.000000 i,ii,j,ek,tt=3 5 6 -2892.901367 0.000000 i,ii,j,ek,tt=3 6 6 21238.716797 1.000000 nodal forces from applied loadsode x-load y-load moment1 0.000000 -192.000000 -204.8000032 0.000000 -192.000000 204.8000033 0.000000 0.000000 0.0000004 0.000000 0.000000 0.000000 i,j,gk1 1 1 241463.421875i,j,gk1 1 2 0.000000i,j,gk1 1 3 7564.800293i,j,gk1 1 4 -238437.500000i,j,gk1 1 5 0.000000i,j,gk1 1 6 0.000000i,j,gk1 1 7 -3025.919922i,j,gk1 1 8 0.000000i,j,gk1 1 9 7564.800293i,j,gk1 1 10 0.000000i,j,gk1 1 11 0.000000i,j,gk1 1 12 0.000000i,j,gk1 2 1 0.000000i,j,gk1 2 2 306642.875000i,j,gk1 2 3 4617.187500i,j,gk1 2 4 0.000000i,j,gk1 2 5 -1442.871094i,j,gk1 2 6 4617.187500i,j,gk1 2 7 0.000000i,j,gk1 2 8 -305200.000000i,j,gk1 2 9 0.000000i,j,gk1 2 10 0.000000i,j,gk1 2 11 0.000000i,j,gk1 2 12 0.000000i,j,gk1 3 1 7564.800293i,j,gk1 3 2 4617.187500i,j,gk1 3 3 44916.000000i,j,gk1 3 4 0.000000i,j,gk1 3 5 -4617.187500i,j,gk1 3 6 9850.000000i,j,gk1 3 7 -7564.800293i,j,gk1 3 8 0.000000i,j,gk1 3 9 12608.000977i,j,gk1 3 10 0.000000i,j,gk1 3 11 0.000000i,j,gk1 3 12 0.000000i,j,gk1 4 1 -238437.500000i,j,gk1 4 2 0.000000i,j,gk1 4 3 0.000000i,j,gk1 4 4 314416.750000i,j,gk1 4 5 -115892.476563i,j,gk1 4 6 4520.158203i,j,gk1 4 8 0.000000i,j,gk1 4 9 0.000000i,j,gk1 4 10 -75979.257813 i,j,gk1 4 11 115892.476563 i,j,gk1 4 12 4520.158203 i,j,gk1 5 1 0.000000i,j,gk1 5 2 -1442.871094 i,j,gk1 5 3 -4617.187500 i,j,gk1 5 4 -115892.476563 i,j,gk1 5 5 184332.921875 i,j,gk1 5 6 -1724.286133 i,j,gk1 5 7 0.000000i,j,gk1 5 8 0.000000i,j,gk1 5 9 0.000000i,j,gk1 5 10 115892.476563 i,j,gk1 5 11 -182890.046875 i,j,gk1 5 12 2892.901367 i,j,gk1 6 1 0.000000i,j,gk1 6 2 4617.187500i,j,gk1 6 3 9850.000000i,j,gk1 6 4 4520.158203i,j,gk1 6 5 -1724.286133 i,j,gk1 6 6 40938.718750 i,j,gk1 6 7 0.000000i,j,gk1 6 8 0.000000i,j,gk1 6 9 0.000000i,j,gk1 6 10 -4520.158203 i,j,gk1 6 11 -2892.901367 i,j,gk1 6 12 10619.358398 i,j,gk1 7 1 -3025.919922 i,j,gk1 7 2 0.000000i,j,gk1 7 3 -7564.800293 i,j,gk1 7 4 0.000000i,j,gk1 7 5 0.000000i,j,gk1 7 6 0.000000i,j,gk1 7 7 3025.919922i,j,gk1 7 8 0.000000i,j,gk1 7 9 -7564.800293 i,j,gk1 7 10 0.000000i,j,gk1 7 11 0.000000i,j,gk1 7 12 0.000000i,j,gk1 8 1 0.000000i,j,gk1 8 2 -305200.000000i,j,gk1 8 4 0.000000i,j,gk1 8 5 0.000000i,j,gk1 8 6 0.000000i,j,gk1 8 7 0.000000i,j,gk1 8 8 305200.000000 i,j,gk1 8 9 0.000000i,j,gk1 8 10 0.000000i,j,gk1 8 11 0.000000i,j,gk1 8 12 0.000000i,j,gk1 9 1 7564.800293i,j,gk1 9 2 0.000000i,j,gk1 9 3 12608.000977i,j,gk1 9 4 0.000000i,j,gk1 9 5 0.000000i,j,gk1 9 6 0.000000i,j,gk1 9 7 -7564.800293i,j,gk1 9 8 0.000000i,j,gk1 9 9 25216.001953i,j,gk1 9 10 0.000000i,j,gk1 9 11 0.000000i,j,gk1 9 12 0.000000i,j,gk1 10 1 0.000000i,j,gk1 10 2 0.000000i,j,gk1 10 3 0.000000i,j,gk1 10 4 -75979.257813 i,j,gk1 10 5 115892.476563 i,j,gk1 10 6 -4520.158203 i,j,gk1 10 7 0.000000i,j,gk1 10 8 0.000000i,j,gk1 10 9 0.000000i,j,gk1 10 10 75979.257813 i,j,gk1 10 11 -115892.476563 i,j,gk1 10 12 -4520.158203 i,j,gk1 11 1 0.000000i,j,gk1 11 2 0.000000i,j,gk1 11 3 0.000000i,j,gk1 11 4 115892.476563 i,j,gk1 11 5 -182890.046875 i,j,gk1 11 6 -2892.901367 i,j,gk1 11 7 0.000000i,j,gk1 11 8 0.000000i,j,gk1 11 9 0.000000i,j,gk1 11 10 -115892.476563i,j,gk1 11 11 182890.046875i,j,gk1 11 12 -2892.901367i,j,gk1 12 1 0.000000i,j,gk1 12 2 0.000000i,j,gk1 12 3 0.000000i,j,gk1 12 4 4520.158203i,j,gk1 12 5 2892.901367i,j,gk1 12 6 10619.358398i,j,gk1 12 7 0.000000i,j,gk1 12 8 0.000000i,j,gk1 12 9 0.000000i,j,gk1 12 10 -4520.158203i,j,gk1 12 11 -2892.901367i,j,gk1 12 12 21238.716797node x-disp y-disp thita1 -0.020768 -0.000749 -0.0041792 -0.021164 -0.014385 0.0078233 -0.000000 -0.000000 0.0000004 0.000000 -0.000000 0.000000reaction nodal forces from the equationsnode x-load y-load moment1 0.249791 -191.991058 -204.7498172 0.253289 -191.827469 204.7061163 94.456879 228.500870 -209.7956244 -94.456657 155.499146 -54.196941element axi-f shear-q moment-m1 94.456924 228.500854 262.488831 -94.456924 155.499146 -28.8833622 228.500870 -94.456879 -209.795624 -228.500870 94.456879 -262.4888313 181.889847 -4.264181 28.883371 -181.889847 4.264181 -54.196941。
有限元方法编程
有限元方法编程简介有限元方法（Finite Element Method, FEM）是一种数值分析方法，用于求解连续介质力学问题。
它将复杂的物理问题离散化为有限数量的简化子问题，并通过构建逼近函数来近似求解。
有限元方法广泛应用于结构力学、流体力学、电磁场等领域。
在本文中，我们将探讨有限元方法的编程实现。
我们将介绍有限元模型的建立、离散化、边界条件的处理以及求解过程等关键步骤。
有限元模型的建立在使用有限元方法求解物理问题之前，首先需要建立一个合适的有限元模型。
这包括选择合适的几何形状和材料属性，并对其进行离散化。
几何形状几何形状是指待求解区域的外部形状，可以是简单几何形状如矩形、圆形，也可以是复杂几何形状如曲线、曲面。
在编程中，我们可以使用各种数学函数或CAD软件来表示几何形状，并将其转换为计算机可识别的格式。
材料属性材料属性是指待求解区域内物质的力学性质，如弹性模量、泊松比等。
在编程中，我们需要将这些属性赋予模型，并将其储存在合适的数据结构中。
离散化离散化是将连续的物理问题转化为离散的子问题。
在有限元方法中，我们通常使用三角形或四边形网格来离散化几何形状。
这些网格被称为有限元网格，每个单元代表一个子问题。
有限元模型的编程实现有限元模型的编程实现主要包括以下几个关键步骤：网格生成、单元属性赋值、边界条件处理以及求解过程。
网格生成在编程中，我们可以使用各种算法来生成有限元网格。
其中最常用的算法是Delaunay三角剖分算法和四边形剖分算法。
这些算法可以根据几何形状和所需精度生成合适的网格，并将其储存在数据结构中。
单元属性赋值每个有限元单元都具有一组属性，包括几何信息和材料信息。
在编程中，我们需要将这些属性赋予每个单元，并将其储存在合适的数据结构中。
这些数据结构通常是数组或矩阵。
边界条件处理边界条件是指在求解过程中所施加的约束条件。
它们可以是位移边界条件、力边界条件或温度边界条件等。
在编程中，我们需要将这些边界条件应用于有限元模型，并将其储存在合适的数据结构中。
有限元编程的c++实现算例
有限元编程的c++实现算例1. #include<>2. #include<>3.4.5. #define ne 3 #define nj 4 #define nz6 #define npj 0 #define npf1 #define nj3 12 #define dd6 #define e0 #definea0 #define i0 #define pi16.17.18. int jm[ne+1][3]={{0,0,0},{0,1,2},{0,2,3},{0,4,3}}; /*gghjghg*/19. double gc[ne+1]={,,,};20. double gj[ne+1]={,,,};21. double mj[ne+1]={,a0,a0,a0};22. double gx[ne+1]={,i0,i0,i0};23. int zc[nz+1]={0,1,2,3,10,11,12};24. double pj[npj+1][3]={{,,}};25. double pf[npf+1][5]={{0,0,0,0,0},{0,-20,,,}};26. double kz[nj3+1][dd+1],p[nj3+1];27. double pe[7],f[7],f0[7],t[7][7];28. double ke[7][7],kd[7][7];29.30.31.36. void jdugd(int);38. void zb(int);39. void gdnl(int);40. void dugd(int);41.42.43. void main()46. int i,j,k,e,dh,h,ii,jj,hz,al,bl,m,l,dl,zl,z,j0;47. double cl,wy[7];48. int im,in,jn;49.50.54. if(npj>0)55. {56. for(i=1;i<=npj;i++)57. { j=pj[i][2];59. p[j]=pj[i][1];60. }61. }62. if(npf>0)63. {64. for(i=1;i<=npf;i++)65. { hz=i;67. gdnl(hz);68. e=(int)pf[hz][3];69. zb(e); for(j=1;j<=6;j++) {72. pe[j]=;73. for(k=1;k<=6;k++) {75. pe[j]=pe[j]-t[k][j]*f0[k];76. }77. }78. al=jm[e][1];79. bl=jm[e][2];80. p[3*al-2]=p[3*al-2]+pe[1]; p[3*al-1]=p[3*al-1]+pe[2];82. p[3*al]=p[3*al]+pe[3];83. p[3*bl-2]=p[3*bl-2]+pe[4];84. p[3*bl-1]=p[3*bl-1]+pe[5];85. p[3*bl]=p[3*bl]+pe[6];86. }87. }89.90. for(e=1;e<=ne;e++) {94. dugd(e); for(i=1;i<=2;i++) {97. for(ii=1;ii<=3;ii++)98. {99. h=3*(i-1)+ii; dh=3*(jm[e][i]-1)+ii; for(j=1;j<=2 ;j++)102. {103. for(jj=1;jj<=3;jj++) {105. l=3*(j-1)+jj; zl=3*(jm[e][j]-1)+jj; dl=zl-dh+ 1; if(dl>0)109. kz[dh][dl]=kz[dh][dl]+ke[h][l]; }111. }112. }113. }114. }115.116. for(i=1;i<=nz;i++) {119. z=zc[i]; kz[z][l]=; for(j=2;j<=dd;j++)122. {123. kz[z][j]=; }125. if((z!=1))126. {127. if(z>dd)128. j0=dd;129. else if(z<=dd)130. j0=z; for(j=2;j<=j0;j++)132. kz[z-j+1][j]=;133. }134. p[z]=; }136.137.138.140. for(k=1;k<=nj3-1;k++)141. {142. if(nj3>k+dd-1) im=k+dd-1;144. else if(nj3<=k+dd-1)145. im=nj3;146. in=k+1;147. for(i=in;i<=im;i++)148. {149. l=i-k+1;150. cl=kz[k][l]/kz[k][1]; jn=dd-l+1;152. for(j=1;j<=jn;j++)153. {154. m=j+i-k;155. kz[i][j]=kz[i][j]-cl*kz[k][m];156. }157. p[i]=p[i]-cl*p[k]; }159. }160.161.162.163.164. p[nj3]=p[nj3]/kz[nj3][1]; for(i=nj3-1;i>=1;i--)166. {167. if(dd>nj3-i+1)168. j0=nj3-i+1;169. else j0=dd; for(j=2;j<=j0;j++)171. {172. h=j+i-1;173. p[i]=p[i]-kz[i][j]*p[h];174. }175. p[i]=p[i]/kz[i][1]; }177. printf("\n");178. printf("_____________________________________________________________\n");179. printf("NJ U V CETA \n"); for(i=1;i<=nj;i++)181. {182. printf(" %-5d % % %\n",i,p[3*i-2],p[3*i-1],p[3*i]);183. }184. printf("_____________________________________________________________\n");185. printf("E N Q M \n");187. for(e=1;e<=ne;e++) {190. jdugd(e); zb(e); for(i=1;i<=2;i++) 193. {194. for(ii=1;ii<=3;ii++)195. {196. h=3*(i-1)+ii;197. dh=3*(jm[e][i]-1)+ii; wy[h]=p[dh];199. }200. }201. for(i=1;i<=6;i++)202. {203. f[i]=;204. for(j=1;j<=6;j++)205. {206. for(k=1;k<=6;k++) {208. f[i]=f[i]+kd[i][j]*t[j][k]*wy[k];209. }210. }211. }212. if(npf>0)213. {214. for(i=1;i<=npf;i++) if(pf[i][3]==e) {217. hz=i;218. gdnl(hz); for(j=1;j<=6;j++) {221. f[j]=f[j]+f0[j];222. }223. }224. }225. printf("%-3d(A) % % %\n",e,f[1],f[2],f[3]); printf(" (B) % % %\n",f[4],f[5],f[6]); }228. return;229. }230.232.236. void gdnl(int hz)237. {238. int ind,e;239. double g,c,l0,d;240.241.242. g=pf[hz][1]; c=pf[hz][2]; e=(int)pf[hz][3];ind=(int)pf[hz][4]; l0=gc[e]; d=l0-c;248. if(ind==1)249. {250. f0[1]=;251. f0[2]=-(g*c*(2-2*c*c/(l0*l0)+(c*c*c)/(l0*l0*l0)))/2; f0[3]=-(g*c*c)*(6-8*c/l0+3*c*c/(l0 *l0))/12;253. f0[4]=;254. f0[5]=-g*c-f0[2];255. f0[6]=(g*c*c*c)*(4-3*c/l0)/(12*l0);256. }257. else258. {259. if(ind==2) {261. f0[1]=;262. f0[2]=(-(g*d*d)*(l0+2*c))/(l0*l0*l0);263. f0[3]=-(g*c*d*d)/(l0*l0);264. f0[4]=;265. f0[5]=(-g*c*c*(l0+2*d))/(l0*l0*l0);266. f0[6]=(g*c*c*d)/(l0*l0);267. }268. else270. f0[1]=-(g*d/l0); f0[2]=;272. f0[3]=;273. f0[4]=-g*c/l0;274. f0[5]=;275. f0[6]=;276. }277. }278. }279.280. void zb(int e)284. {285. double ceta,co,si;286. int i,j;287. ceta=(gj[e]*pi)/180; co=cos(ceta); 289. si=sin(ceta);290. t[1][1]=co; t[1][2]=si;292. t[2][1]=-si;293. t[2][2]=co;294. t[3][3]=;295. for(i=1;i<=3;i++)296. {297. for(j=1;j<=3;j++) {299. t[i+3][j+3]=t[i][j];300. }301. }302. }303.304.305.306. void jdugd(int e)310. {311. double A0,l0,j0;312. int i;314.315.316. A0=mj[e]; l0=gc[e]; j0=gx[e];320.321. for(i=0;i<=6;i++)322. for(j=0;j<=6;j++) kd[i][j]=;324.325. kd[1][1]=e0*A0/l0;326. kd[2][2]=12*e0*j0/pow(l0,3);327. kd[3][2]=6*e0*j0/pow(l0,2);328. kd[3][3]=4*e0*j0/l0;329. kd[4][1]=-kd[1][1];330. kd[4][4]=kd[1][1];331. kd[5][2]=-kd[2][2]; kd[5][3]=-kd[3][2];333. kd[5][5]=kd[2][2];334. kd[6][2]=kd[3][2];335. kd[6][3]=2*e0*j0/l0;336. kd[6][5]=-kd[3][2];337. kd[6][6]=kd[3][3];338.339. for(i=1;i<=6;i++)340. for(j=1;j<=i;j++) kd[j][i]=kd[i][j];342. }343.344.345. void dugd(int e)349. {350. int i,k,j,m;351. jdugd(e); zb(e); for(i=1;i<=6;i++)354. {355. for(j=1;j<=6;j++)356. {357. ke[i][j]=;358. for(k=1;k<=6;k++)359. {360. for(m=1;m<=6;m++)361. {362. ke[i][j]=ke[i][j]+t[k][i]*kd[k][m]*t[m][j]; } 364. }365. }366. }367. }368.369.370.372.373.374. </></>。
平面桁架有限元C#编程
1题目结构如图所示： 杆的弹性模量E 为200000Mpa ，横截面面积A 为3250mm 2。
图 1 桁架示意图2实验材料PC 机一台，Microsoft Visual Studio 软件，Ansys 软件。
3实验原理（1）桁架结构特点桁架结构中的桁架指的是桁架梁，是格一种梁式结构。
桁架结构常用于大跨度的厂房、展览馆、体育馆和桥梁等公共建筑中。
由于大多用于建筑的屋盖结构，桁架通常也被称作屋架。
结构上由光滑铰链连接，载荷只作用于节点处，只约束线位移，各杆只有轴向拉伸和压缩。
（2）平面桁架有限元分析1、单元分析局部坐标系中的干单元如图所示：图 2 局部坐标系中的杆单元以下公式描述了整体位移和局部位移之间的关系：U=Tu 其中U=[ U ix U iy U jx U jy ]，T=[cos θ−sin θ00sin θcos θ0000cos θ−sin θ00sin θcos θ],u=[u ix u iy u jx u jy ]U 和u 分别代表整体坐标系和局部坐标系XY 系和局部坐标系xy 下节点i 和节点j 的位移。
T 是变形从局部坐标转换到整体坐标系下的变换阵，类似的局部力和整体力也有以下关系：F=Tf其中F=[ F ixF iy F jx F jy ] ，是整体坐标系下施加在节点i 和j 上的力的分量而且f=[ f ix f iy f jx f jy ]，代表局部坐标系下施加在节点i和j上的分量。
在假设的二力杆条件下，杆只能沿着局部坐标系的x方向变形，内力也总是沿着局部坐标系x的方向，因此将y方向的位移设置为0，局部坐标系下内力和位移通过刚度矩阵有如下关系：[f ixf iyf jxf jy]=|k0−k00000−k0k00000|=[U ixU iyU jxU jy]这里k=k eq=AE/L，写成矩阵形式有：f=Ku将f和u替换成F和U有：T-1F=KT-1U将方程两边乘以T得到：F=TKT-1U其中T-1是变换矩阵T的逆矩阵，替换方程中的TKT-1和U矩阵的值，相乘后得到：[F ixF iy F jx F jy]= k[cos2θsinθcosθ−cos2θ−sinθcosθsinθcosθsin2θ−sinθcosθ−sin2θ−cos2θ−sinθcosθcos2θsinθcosθ−sinθcosθ−sin2θsinθcosθsin2θ][U ixU iyU jxU jy]上述方程代表了施加外力、单元刚度矩阵和任意单元节点的整体位移之间的关系。
有限单元法的编程实现，内附Fortran源代码
作业2：有限单元法的编程实现2.1有限元概述有限元分析的基本概念是将求解域离散为若干子域，并通过它们边界上的节点相互联结成为组合体，对每一单元假定一个合适的近似解，然后推导求解这个域总的满足条件，从而得到问题的解。
这个解只是近似解，因为实际问题被较简单的问题所代替。
由于大多数实际问题难以得到准确解，而有限元不仅计算精度高，而且能适应各种复杂形状，因而成为行之有效的工程分析手段。
对于不同物理性质和数学模型的问题，有限元求解法的基本步骤是相同的，只是具体公式推导和运算求解不同。
有限元求解问题的基本步骤通常为：2.1.1确定求解域及求解域的离散化，即子域的剖分根据实际问题近似确定求解域的物理性质和几何区域。
将求解域近似为具有不同有限大小和形状且彼此相连的有限个单元组成的离散域，习惯上称为有限元网络划分，求解域的离散化是有限元法的核心技术之一。
选择单元的形式，确定单元数，节点数及单元及节点的编号。
2.1.2单元分析一个具体的物理问题通常可以用一组包含问题状态变量边界条件的微分方程式表示，为适合有限元求解，通常将微分方程化为等价的泛函形式。
作为弹性力学微分方程的等效积分的形式，虚位移原理与虚应力原理分别是平衡方程与力的边界条件和几何方程与位移边界条件的等效积分形式。
在导出它们的的过程中都未涉及到物理方程所以它们适用于线弹性、非线性线弹性及弹塑性的问题。
现有的有限元计算多采用以位移为未知量的形式，可用虚位移原理来描述其平衡方程，其矩阵形式为：(u )u 0T T T V S f dV TdS σδεσδδ--=⎰⎰ (2.1) 我们需要对单元构造一个适合的近似解，即推导有限单元的格式，其中包括选择合理的单元坐标系，建立单元试函数，以某种方法给出单元各状态变量的离散关系，从而形成单元矩阵（结构力学中称刚度阵或柔度阵）。
这里以平面问题为例，单元位移与节点位移的关系表示为：[,,...]e i e e i i i j j a u u N a N N a Na ⎧⎫⎪⎪====⎨⎬⎪⎪⎩⎭∑ (2.2) 应变与节点位移的关系为：[][]e e e e i j i j Lu LNa L N N a B B a Ba ε===== (2.3)应力与节点位移的关系为：e e D DBa Sa σε=== (2.4)单元刚度矩阵的表达式为：T ee K B DBtdxdy Ω=⎰ (2.5) 单元等效节点载荷矩阵表示为T T e ee e ef s P P P N ftdxdy N TtdS ΩΩ=+=+⎰⎰ (2.6) 2.1.3 整体分析将单元总装形成离散域的总矩阵方程（联合方程组），反映对近似求解域的离散域的要求，即单元函数的连续性要满足一定的连续条件。
有限元C++编程实践
基于有限元算法的编程实践学号：2011043010031姓名：廖校毅电子科技大学物理电子学院【摘要】本文就有限元算法在电磁场分析中的应用展开研究与实践，从电磁场的Maxwell方程出发,根据电磁场的边值问题及变分公式建立了有限元方程组。
具体在实践中，将这些知识运用到C++语言和Matlab中，并将这两种语言有机结合，编程并实现二维FEM。
程序最后通过计算含两种介质的电位槽电位分布来验证其正确性。
关键词: 有限元变分法C++ MatlabThe Programming Practice Based on The Finite Element Algorithm Student ID：2011043010031Name：Liao Xiaoyi University of Electronic Science and technology &School of PhysicalElectronicsAbstract In this paper, we take the application of finite element method in electromagnetic field analysis into research and practice. Starting from Maxwell equations of electromagnetic field, the electromagnetic field boundary value problems and variational formula established the finite element equations. In specific practice, this knowledge will be applied to the C++ language and Matlab, and the organic combination of two languages, programming and implementation of two-dimensional FEM. Finally, through the program to verify the validity of the calculation of potential distribution in channel potential containing two kinds of medium.Key words The finite element method The variational method C++ Matlab一、前言在数学中，有限元法（FEM，Finite Element Method）是一种为求得偏微分方程边值问题近似解的数值技术。