带裂纹的动态有限元程序（水平裂纹在中间高度处，从左到中点--即半无限长）
来源：互联网 发布：软件委托开发合同 编辑：程序博客网 时间：2024/04/30 10:23
Crack_DynamicFEM.java
package scu.edu.fem;
class Crack_DynamicETNode{ //单元结点结构体
double x=0,y=0; //单元结点坐标
int number=0; //单元结点在总体区域划分中的编号
}
class Crack_DynamicElementTriangle{ //三角形单元结构体
Crack_DynamicETNode nd[]=new Crack_DynamicETNode[]{new Crack_DynamicETNode(),new Crack_DynamicETNode(),new Crack_DynamicETNode()}; //存储相对应的总体结点号
double a[]=new double[3];
double b[]=new double[3];
double c[]=new double[3]; //基函数的系数
double A; //单元面积
// double Aij[3][3]; //单元矩阵
double sigma11,sigma22,sigma33,sigma12; //正应力，及剪切应力
double stain11,strain22,strain12; //正应变，及剪切应变
//jimmy 改成6*6的二维
double Aij[][]= new double[6][6]; //单元矩阵(K)
double M[][]=new double[6][6];; //单元矩阵(M)
double C[][]=new double[6][6];; //单元矩阵(C)
double fi[]=new double[3];; //单元上的积分值
}
public class Crack_DynamicFEM {
// int i; //单元的循环指标
// int j;
// int k;
//int id;
static double dBig=1000000000000000000000000000000f;
//jimmy change for test
//const int nx=1,ny=1; //x,y方向划分网格数,三角形单元个数=nx*ny*2
//jimmy change for test 必须为奇数
final static int nx=5,ny=5; //x,y方向划分网格数,三角形单元个数=nx*ny*2,这里对裂纹假设是5*5
//const int nx=21,ny=21; //x,y方向划分网格数,三角形单元个数=nx*ny*2
final static double Lx=10; //矩形区域的两边长
final static double Ly=10;
static double D; //单元矩阵行列式值
static double E=10000000f; //模量
static double t=1.0; //厚度
static double u=0.3; //泊松比
static double rou=1.0; //密度
static double Y=0.00000; //阻尼系.0
static double T=1;
static double dataT=0.01; //时间步长
static double gama=0.5; //参数gama>=0.5
static double bata=1;//0.25*(0.5+gama)*(0.5+gama) //参数bata
static double c0=1/(bata*dataT*dataT); //参数c0
static double c1=gama/(bata*dataT); //参数c1
static double c2=1/(bata*dataT); //参数c2
static double c3=1/(bata*2)-1; //参数c3
static double c4=gama/bata-1; //参数c4
static double c5=dataT/2*(c4-1); //参数c5
static double c6=dataT*(1-gama); //参数c6
static double c7=gama*dataT; //参数c7
static double pMatrix[]; //总体K矩阵指针
static double totalM[]; //总体质量M矩阵指针
static double totalC[]; //总体阻尼C矩阵指针
static double dalta[]; //位移
static double dalta1[]; //速度
static double dalta2[]; //加速度
static double pMf[]; //f向量指针
static double tmp[];
static Crack_DynamicElementTriangle[] pE; //单元三角形结构体数组指针
double ai,bi,ci; //基函数的系数
//-----------------------------------------------------
double fn(int n , double x[])//被积函数,高斯积分函数的参数
{
return(ai+bi*x[0]+ci*x[1]);
}
//---------- 全选主元高斯消去法 ------------------------------
// a 体积为n*n的双精度实型二维数组，方程组系数矩阵，返回时将被破坏
// b 长度为n的双精度实型一维数组，方程组右端的常数向量，返回方程组的解向量
// n 整型变量，方程组的阶数
//--------------------------------------------------------------
static int Gauss(double a[],double b[],int n)
{
int []js;
int l,k,i,j,is=0,p,q;
double d,t;
js=new int[n];
l=1;
for(k=0;k<=n-2;k++)
{
d=0.0;
for(i=k;i<=n-1;i++)
for(j=k;j<=n-1;j++)
{
t=Math.abs(a[i*n+j]);
if(t>d) { d=t; js[k]=j; is=i;}
}
if(d+1.0==1.0) l=0;
else
{
if(js[k]!=k)
for(i=0;i<=n-1;i++)
{
p=i*n+k; q=i*n+js[k];
t=a[p]; a[p]=a[q]; a[q]=t;
}
if(is!=k)
{
for(j=k;j<=n-1;j++)
{
p=k*n+j; q=is*n+j;
t=a[p]; a[p]=a[q]; a[q]=t;
}
t=b[k]; b[k]=b[is]; b[is]=t;
}
}
if(l==0)
{
System.out.println("Gauss funtion failed 1...\n");
return(0);
}
d=a[k*n+k];
for(j=k+1;j<=n-1;j++)
{
p=k*n+j; a[p]=a[p]/d;
}
b[k]=b[k]/d;
for(i=k+1;i<=n-1;i++)
{
for(j=k+1;j<=n-1;j++)
{
p=i*n+j;
a[p]=a[p]-a[i*n+k]*a[k*n+j];
}
b[i]=b[i]-a[i*n+k]*b[k];
}
}
d=a[(n-1)*n+n-1];
if(Math.abs(d)+1.0==1.0)
{
//free(js);
System.out.println("Gauss funtion failed 2...\n");
return(0);
}
b[n-1]=b[n-1]/d;
for(i=n-2;i>=0;i--)
{
t=0.0;
for(j=i+1;j<=n-1;j++)
t=t+a[i*n+j]*b[j];
b[i]=b[i]-t;
}
js[n-1]=n-1;
for(k=n-1;k>=0;k--)
if(js[k]!=k)
{
t=b[k]; b[k]=b[js[k]]; b[js[k]]=t;
}
return 1;
}
public static void main(String[] args) {
//-------------- 全局变量 ---------------------------
Crack_DynamicFEM staticFEM=new Crack_DynamicFEM();
final int iNode=(nx+1)*(nx+1); //结点个数
//为总体矩阵,三角形单元数组,f函数向量分配存储内存
pMatrix=new double[2*iNode*2*iNode];
totalM=new double[2*iNode*2*iNode];
totalC=new double[2*iNode*2*iNode];
pE=new Crack_DynamicElementTriangle[nx*ny*2];
for(int i=0;i<nx*ny*2;i++)
pE[i]=new Crack_DynamicElementTriangle();
//jimmy 修改为2*iNode，因为有x,y两个方向载荷或约束
pMf=new double[2*iNode];
tmp=new double[2*iNode];
dalta=new double[2*iNode];
dalta1=new double[2*iNode];
dalta2=new double[2*iNode];
//初始化值为0,因为下面要累加总体矩阵
for(int i=0;i<2*iNode*2*iNode;i++)
totalC[i]=totalM[i]=pMatrix[i]=0.0;
for(int i=0;i<2*iNode;i++)
dalta2[i]=dalta[i]=dalta1[i]=pMf[i]=0.0;
try{
//------ 计算得到网格的信息 -----------
for(int j=0;j<nx;j++)
for(int i=0;i<ny;i++)
{
//for the first triangle in the rectangle
pE[i*2+j*ny*2].nd[0].x=(Lx/nx)*i;
pE[i*2+j*ny*2].nd[0].y=(Ly/ny)*j;
pE[i*2+j*ny*2].nd[0].number=i+j*(nx+1); //NO.0
pE[i*2+j*ny*2].nd[1].x=(Lx/nx)*(i+1);
pE[i*2+j*ny*2].nd[1].y=(Ly/ny)*(j+1);
pE[i*2+j*ny*2].nd[1].number=i+1+(nx+1)*(j+1); //NO.1
//加边缘裂纹（当划分为5*5时）在单元i==0&&j==2 ，及i==0, j==3两单元加裂纹
if(i==0&&j==2){
pE[i*2+j*ny*2].nd[2].x=(Lx/nx)*i;
pE[i*2+j*ny*2].nd[2].y=(Ly/ny)*(j+1)-0.1;
pE[i*2+j*ny*2].nd[2].number=i+(nx+1)*(j+1); //NO.2
}
else{
pE[i*2+j*ny*2].nd[2].x=(Lx/nx)*i;
pE[i*2+j*ny*2].nd[2].y=(Ly/ny)*(j+1);
pE[i*2+j*ny*2].nd[2].number=i+(nx+1)*(j+1); //NO.2
}
//for the second triangle in the rectangle
if(i==0&&j==3){
pE[i*2+j*ny*2+1].nd[0].x=(Lx/nx)*i;
pE[i*2+j*ny*2+1].nd[0].y=(Ly/ny)*j+0.1;
pE[i*2+j*ny*2+1].nd[0].number=i+j*(nx+1); //NO.0
}
else{
pE[i*2+j*ny*2+1].nd[0].x=(Lx/nx)*i;
pE[i*2+j*ny*2+1].nd[0].y=(Ly/ny)*j;
pE[i*2+j*ny*2+1].nd[0].number=i+j*(nx+1); //NO.0
}
pE[i*2+j*ny*2+1].nd[1].x=(Lx/nx)*(i+1);
pE[i*2+j*ny*2+1].nd[1].y=(Ly/ny)*j;
pE[i*2+j*ny*2+1].nd[1].number=i+j*(nx+1)+1; //NO.1
pE[i*2+j*ny*2+1].nd[2].x=(Lx/nx)*(i+1);
pE[i*2+j*ny*2+1].nd[2].y=(Ly/ny)*(j+1);
pE[i*2+j*ny*2+1].nd[2].number=i+1+(nx+1)*(j+1); //NO.2
}
System.out.println("计算基函数系数值...\n");
int j=1,k=2;
for(int id=0;id<nx*ny*2;id++)
{
for( int 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;
}
}System.out.println("OK!\n");
System.out.println("计算单元有限元特征式系数矩阵...\n");
//--------------------------------------计算单元有限元特征式系数矩阵可以不再分两个三角形循环
//Aij即是Et/[4*(1-u*u)*A]*一个2*2Matrix
int l,m;
for(int i=0;i<nx*ny*2;i++)
{ for(l=0;l<3;l++)
for(m=0;m<3;m++) // Respaired
{
pE[i].Aij[2*l][2*m]=( pE[i].b[l]*pE[i].b[m]+ pE[i].c[l]*pE[i].c[m]*(1-u)/2 ) *E*t/(2*(1-u*u));
pE[i].Aij[2*l][2*m+1]=( pE[i].b[l]*pE[i].c[m]*u + pE[i].c[l]*pE[i].b[m]*(1-u)/2)*E*t/(2*(1-u*u));
pE[i].Aij[2*l+1][2*m]=( pE[i].c[l]*pE[i].b[m]*u + pE[i].b[l]*pE[i].c[m]*(1-u)/2 ) *E*t/(2*(1-u*u));
pE[i].Aij[2*l+1][2*m+1]=( pE[i].c[l]*pE[i].c[m]*u + pE[i].b[l]*pE[i].b[m]*(1-u)/2 ) *E*t/(2*(1-u*u));
}
System.out.println("\n");
// System.out.println("%10.2f\t",pE[i].Aij[2*l+1][2*m]);
// System.out.println("%10.2f\t",pE[i].Aij[2*l+1][2*m+1]);
}
for(int i=0;i<nx*ny*2;i++)
{ for(l=0;l<3;l++)
for(m=0;m<3;m++) // Respaired
{
if(l==m)
{
pE[i].M[2*l][2*m]=(rou)*t*pE[i].A/6;
pE[i].M[2*l][2*m+1]=0.0;
pE[i].M[2*l+1][2*m]=0.0;
pE[i].M[2*l+1][2*m+1]=(rou)*t*pE[i].A/6;
pE[i].C[2*l][2*m]=Y*t*pE[i].A/6;
pE[i].C[2*l][2*m+1]=0.0;
pE[i].C[2*l+1][2*m]=0.0;
pE[i].C[2*l+1][2*m+1]=Y*t*pE[i].A/6;
}else
{
pE[i].M[2*l][2*m]=(rou)*t*pE[i].A/12;
pE[i].M[2*l][2*m+1]=0.0;
pE[i].M[2*l+1][2*m]=0.0;
pE[i].M[2*l+1][2*m+1]=(rou)*t*pE[i].A/12;
pE[i].C[2*l][2*m]=Y*t*pE[i].A/12;
pE[i].C[2*l][2*m+1]=0;
pE[i].C[2*l+1][2*m]=0;
pE[i].C[2*l+1][2*m+1]=Y*t*pE[i].A/12;
}
}
System.out.println("\n");
// System.out.println("%10.2f\t",pE[i].Aij[2*l+1][2*m]);
// System.out.println("%10.2f\t",pE[i].Aij[2*l+1][2*m+1]);
}
System.out.println("打印单元特征式系数矩阵...\n");
//fSystem.out.println(wfp,"计算得各结点上的温度值为:\n");
for(int i=0;i<nx*ny*2;i++)
{
System.out.println("计算得单元"+i+"上的K矩阵为:\n");
for(l=0;l<6;l++)
{
for(m=0;m<6;m++)
{
// if(abs(pE[i].Aij[l][m]-0.0)<0.1)
// System.out.println("0.00\t");
// else
System.out.println("\t"+pE[i].Aij[l][m]);
}
System.out.println("\n");
// System.out.println("%6.2f\t",pE[i].Aij[2*l+1][2*m]);
// System.out.println("%6.2f\t",pE[i].Aij[2*l+1][2*m+1]);
}
}
//fSystem.out.println(wfp,"%d %f\n",i+1,pMf[i]);
// System.out.println("写特征式系数矩阵到文件OK!\n");
// fclose(wfp);
System.out.println("计算积分值,填充到f函数向量数组...\n");
//---------------------------计算三角形单元f函数向量不用多重积分公式,直接代用单元三角形面积的三分之一
int idx=0;
for(int i=0;i<2*nx*ny;i++)
{ for(idx=0;idx<3;idx++)
{
pE[i].fi[idx]=(0.5)*(Lx/nx)*(Ly/ny); // Respaired
}
}
// System.out.println("pE[0].fi[0]=%f\n",pE[0].fi[0]);
// System.out.println("pE[0].fi[1]=%f\n",pE[0].fi[1]);
// System.out.println("pE[0].fi[2]=%f\n",pE[0].fi[2]);
// System.out.println("OK!\n");
//单元矩阵元素累加到总体矩阵相应的位置上,共四个节点，两个单元
// 2 |-------|3
// | I . |
// | . II |
// 0 |-------|1
System.out.println("单元矩阵元素累加到总体矩阵相应的位置上...\n");
for(idx=0;idx<nx*ny*2;idx++)
for(int i=0;i<3;i++)
{
for(j=0;j<3;j++)
{ pMatrix[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)] +=pE[idx].Aij[2*i][2*j];
pMatrix[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)+1]+=pE[idx].Aij[2*i][2*j+1];
pMatrix[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)+2*iNode] +=pE[idx].Aij[2*i+1][2*j];
pMatrix[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)+2*iNode+1] +=pE[idx].Aij[2*i+1][2*j+1];
totalM[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)] +=pE[idx].M[2*i][2*j];
totalM[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)+1]+=pE[idx].M[2*i][2*j+1];
totalM[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)+2*iNode] +=pE[idx].M[2*i+1][2*j];
totalM[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)+2*iNode+1] +=pE[idx].M[2*i+1][2*j+1];
totalC[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)] +=pE[idx].C[2*i][2*j];
totalC[2*(pE[idx].nd[i].number*iNode+pE[idx].nd[j].number)+1]+=pE[idx].C[2*i][2*j+1];
totalC[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)+2*iNode] +=pE[idx].C[2*i+1][2*j];
totalC[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)+2*iNode+1] +=pE[idx].C[2*i+1][2*j+1];
//pMatrix[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)] += (totalM[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)]+totalC[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)]);
//pMatrix[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)+1]+=( totalM[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)+1]+totalC[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)+1]);
//pMatrix[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)+2*iNode] += (totalM[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)+2*iNode]+totalC[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)+2*iNode]);
//pMatrix[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)+2*iNode+1] +=(totalM[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)+2*iNode+1]+totalC[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)+2*iNode+1]);
//无阻尼
pMatrix[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)] += (totalM[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)]);
pMatrix[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)+1]+=( totalM[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)+1]);
pMatrix[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)+2*iNode] += (totalM[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)+2*iNode]);
pMatrix[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)+2*iNode+1] +=(totalM[2*(pE[idx].nd[i].number*2*iNode+pE[idx].nd[j].number)+2*iNode+1]);
pMf[pE[idx].nd[i].number ]+=pE[idx].fi[i];
}
}
/*
for(i=0;i<iNode*2;i++)
{
for(j=0;j<iNode*2;j++)
{
System.out.println("%6d\t",pMatrix[i*iNode*2+j]);
}
System.out.println("%6d\t",pMf[i*iNode*2+j]);
System.out.println("\n");
}
*/
// System.out.println("总体矩阵到文件...\n");
//文件操作
// FILE *wfp;
// if((wfp=fopen("gobalMatrix.txt","w"))==NULL)
// System.out.println("Cann't open the file... ");
for(idx=0;idx<nx*ny*2;idx++)
{
System.out.println("\r\n");
for(int i=0;i<3;i++)
{
for(int h=0;h<3;h++)
{
//System.out.println("%6.2f\t"+pE[idx].Aij[2*i][2*h]);
//System.out.println("%6.2f\t"+pE[idx].Aij[2*i][2*h+1]);
}
//System.out.println("\n");
//System.out.println("%6.2f\t"+pE[idx].Aij[2*i+1][2*2]);
//System.out.println("%6.2f\t"+pE[idx].Aij[2*i+1][2*2+1]);
}
}
//fSystem.out.println(wfp,"%d %f\n",i+1,pMf[i]);
//System.out.println("\r\n总体矩阵到文件OK!\n");
//fclose(wfp);
//double dBig=pow(10,20); //边界条件对角线扩大法处理所用的大数
// double Ur=0.0; //边界条件1(边界条件2通过Calerkin弱解表达式自动满足)
double Ur=1000 ; //外力，牛
//jimmy 上下拉伸
//上端2,3节点，Y方向拉力 索引为5,7
for(int i=0;i<nx+1;i++)
{
pMf[(nx+1)*nx*2+2*i+1]=Ur;
}
//下端0,1节点，Y方向拉力为负 ，索引为1,3
// 位移 ，索引
//u0, 0
// v0, 1
/// u1, 2
// v1, 3
// u2, 4
// v2, 5
// u3, 6
// v3, 7
//设置为y方向为0位移，参考王勋成的书
for(j=0;j<nx+1;j++)
{
pMatrix[(2*j)*2*iNode+2*j]*=dBig;
pMatrix[(2*j+1)*2*iNode+2*j+1]*=dBig;
pMf[2*j]=0;
pMf[2*j+1]=0;
}
//右端 for(j=0;j<nx+1;j++)
// { i=(nx+1)*(j+1)-1;
// pMatrix[i*iNode+i]*=dBig;
// pMf[i]=pMatrix[i*iNode+i]*Ur;
//}
System.out.println("调用全选主元高斯消去法函数解方程组...\n");
//if((wfp=fopen("result.txt","w"))==NULL)
// System.out.println("Cann't open the file... ");
//假设一
//初始加速度为确定值
// for(int i=iNode-1;i<iNode;i++)
// dalta2[2*i-1]=-0.1;
//
//假设二（与假设一稍有不同）
//第一步求t=0加速度的，结果在pMf中
//假设t=0，速度和位移为0
for(int i=0;i<2*iNode;i++)
tmp[i]=pMf[i];
Gauss(pMatrix,tmp,2*iNode);
for(int i=0;i<2*iNode;i++)
dalta[i]=tmp[i];
// Gauss(totalM,pMf,2*iNode);
// for(int i=iNode-(nx+1);i<iNode;i++)
// dalta2[2*i-1]=-pMf[i];
//假设三（与假设一,二稍有不同，这里给定初始速度）
//for(int i=(2*iNode-1);i<2*iNode;i++)
//dalta1[i]=-1;
//增加时间步长的内容
for(int step=1;step<T/dataT;step++)
{
//if((wfp=fopen("result2.txt","w+"))==NULL)
//System.out.println("Cann't open the file... ");
for(int i=0;i<2*iNode;i++)
for(j=0;j<2*iNode;j++)
{
//pMf[i]=totalM[i*2*iNode+j]*(c0*dalta[j]+c2*dalta1[j]+c3*dalta2[j])+totalC[i*2*iNode+j]*(c1*dalta[j]+c4*dalta1[j]+c5*dalta2[j]);
pMf[i]=pMf[i]+totalM[i*2*iNode+j]*(c0*dalta[j]+c2*dalta1[j]+c3*dalta2[j]);
}
Gauss(pMatrix,pMf,2*iNode); //调用全选主元高斯消去法函数解方程组
for(int i=0;i<2*iNode;i++)
{
dalta2[i]=c0*(pMf[i]-dalta[i])-c2*dalta1[i]-c3*dalta2[i];
dalta1[i]=dalta1[i]+c6*dalta2[i]+c7*dalta2[i];
dalta[i]=pMf[i];
}
//System.out.println("节点编号："+iNode+"\t"+pMf[2*iNode-2]+"\t"+pMf[2*iNode-1]);
//System.out.println(pMf[2*iNode-2]+"\t"+pMf[2*iNode-1]);
//System.out.println(pMf[2*iNode-1]);
//int i=49;
pE[23].sigma11=( pE[23].b[0]*pMf[26]+ u*pE[23].c[0]*pMf[27]+pE[23].b[1]*pMf[40]+ u*pE[23].c[1]*pMf[41]+pE[23].b[2]*pMf[38]+ u*pE[23].c[2]*pMf[39])*E*t/(2*pE[23].A*(1-u*u));
pE[23].sigma22=( u*pE[23].b[0]*pMf[26]+ pE[23].c[0]*pMf[27]+u*pE[23].b[1]*pMf[40]+ pE[23].c[1]*pMf[41]+u*pE[23].b[2]*pMf[38]+pE[23].c[2]*pMf[39])*E*t/(2*pE[23].A*(1-u*u));
//pE[i].sigma12=( pE[i].c[l]*pE[i].b[m]*u + pE[i].b[l]*pE[i].c[m]*(1-u)/2 ) *E*t/(2*(1-u*u));
System.out.println(pMf[2*iNode-1]);
//System.out.println(pMf[2*iNode-1]+"\t"+pE[23].sigma11+"\t"+pE[23].sigma22);
//System.out.println("\r\n");
// fclose(wfp);
//System.out.println("step="+step);
}
}
catch(Exception e)
{
e.printStackTrace();
System.out.println("Error occured...\n");
}
//释放总体矩阵和三角形单元数组占用内存
//free(pMf); free(pE); free(pMatrix);
System.out.println("Hello World!\n");
return ;
}
}
- 带裂纹的动态有限元程序（水平裂纹在中间高度处，从左到中点--即半无限长）
- 转子或轴裂纹的诊断要点
- 恳戮肯饰豢翰剖谄惨裂纹撤菏菇逗
- PS打造裂纹文字
- 详细分析裂纹止裂试验装置
- 数据处理--文章《按照朱老师修改边框_不对称共线裂纹的应力强度因子的相互影响关系研究》
- distmsh网格划分之长方形含边界裂纹
- 有限元算法-4：有限元代码（由Shell63，Solid45，Fluid30组合在一起形成的耦合程序）
- 有限元方法基础入门教程（一维弹性问题的有限元程序）
- 基于并行EBE-CG方法的有限元求解程序（从我的毕业论文中节选出来的）
- 特征值问题的有限元MATLAB程序（一维）
- 带约束有限元网格的平滑优化
- IOS--如何从第N级界面返回到一级界面（即添加在UITabbarController上的UIController）
- 程序越长水平越高吗？
- 动态有限元程序（二维杆，上端拉伸后放手，下端固定）
- UFT 12的破解（即无限使用）
- 从OpenGL 1.x 到 2.x的迁移（即从固定管线到可编程管线的迁移）
- 输入无限长的字符
- OnCtrlColor,OnDrawItem,DrawItem,OnPaint之间关系
- java WEB中的定时器
- android获取应用签名信息
- FAQ_23 设置 Toast 显示时间
- 使用Zend Studio和火狐调试PHP代码
- 带裂纹的动态有限元程序（水平裂纹在中间高度处，从左到中点--即半无限长）
- Coverity的安装使用及常见错误
- DirectX 11 framebuffer capture (C , no Win32 or D3DX)
- Delphi 报表 Quickreport使用
- IE6、IE7下解决Select框不随滚动条滚动的问题
- 构造函数的执行与virtual
- 笔记1
- 二维动态有限差分程序
- VS2008 可执行文件路径选择
