DynamicFDM.java
package scu.edu.fem;
class FDMETNode_2D{ //单元结点结构体
double x=0,y=0; //单元结点坐标
int number=0; //单元结点在总体区域划分中的编号
double accX; //节点加速度
double accY; //节点加速度
double velocityX; //节点速度矩阵
double velocityY; //节点速度矩阵
double displacementX; //节点速度矩阵
double displacementY; //节点速度矩阵
}
class FDMElementTriangle_2D implements Cloneable{ //四边形单元结构体
FDMETNode_2D nd[]=new FDMETNode_2D[]{new FDMETNode_2D(),new FDMETNode_2D(),new FDMETNode_2D(),new FDMETNode_2D()}; //存储相对应的总体结点号
double a[]=new double[3];
double b[]=new double[3];
double c[]=new double[3]; //基函数的系数
double A; //单元面积
double sigmaX,sigmaY,sigmaZ,sigmaXY; //正应力，及剪切应力
double strainX,strainY,strainXY; //正应变，及剪切应变
double fi[]=new double[3]; //单元上的积分值
double strainRatioX; //单元应变率
double strainRatioY; //单元应变率
double strainRatioXY; //单元应变率XY
double p,Sx,Sy,Sxy;
// clone(){
//
// }
@Override
protected Object clone() throws CloneNotSupportedException {
// TODO Auto-generated method stub
return super.clone();
}
}
public class DynamicFDM_2D {
// 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=100,ny=100; //x,y方向划分网格数,三角形单元个数=nx*ny*2
//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=10000f; //模量
static double t=1.0; //厚度
static double u=0.3; //泊松比
static double rou=1.0; //密度
static double Y=0.1; //阻尼系数
static double T=10;
static double dataT=0.001; //时间步长
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 total_accX[]; //节点加速度指针
static double total_accY[]; //节点加速度指针
static double total_velocityX[]; //节点速度矩阵指针
static double total_velocityY[]; //节点速度矩阵指针
static double total_displacementX[]; //节点速度矩阵指针
static double total_displacementY[]; //节点速度矩阵指针
static double total_strainRatioX[]; //总体应变率矩阵指针
static double total_strainRatioY[]; //总体应变率矩阵指针
//总外力Fx,Fy
static double Fx[]; //f向量指针
static double Fy[]; //f向量指针
static double tmp[];
static FDMElementTriangle_2D[] pE; //单元三角形结构体数组指针
double ai,bi,ci; //基函数的系数
//-----------------------------------------------------
double fn(int n , double x[])//被积函数,高斯积分函数的参数
{
return(ai+bi*x[0]+ci*x[1]);
}
public static void main(String[] args) {
//-------------- 全局变量 ---------------------------
DynamicFDM_2D staticFEM=new DynamicFDM_2D();
final int iNode=(nx+1)*(nx+1); //结点个数
pE=new FDMElementTriangle_2D[nx*ny];
for(int i=0;i<nx*ny;i++)
pE[i]=new FDMElementTriangle_2D();
total_displacementX=new double[iNode];
total_velocityX=new double[iNode];
total_accX=new double[iNode];
total_displacementY=new double[iNode];
total_velocityY=new double[iNode];
total_accY=new double[iNode];
total_strainRatioX=new double[iNode];
total_strainRatioY=new double[iNode];
//jimmy 修改为2*iNode，因为有x,y两个方向载荷或约束
Fx=new double[iNode];
Fy=new double[iNode];
//tmp=new double[2*iNode];
for(int i=0;i<iNode;i++)
Fx[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,与总体网格保持一致
//2-----3
//| |
//0-----1
pE[i+j*ny].nd[0].x=(Lx/nx)*i;
pE[i+j*ny].nd[0].y=(Ly/ny)*j;
pE[i+j*ny].nd[0].number=i+j*(nx+1); //NO.0
pE[i+j*ny].nd[1].x=(Lx/nx)*(i+1);
pE[i+j*ny].nd[1].y=(Ly/ny)*j;
pE[i+j*ny].nd[1].number=i+j*(nx+1)+1; //NO.1
pE[i+j*ny].nd[2].x=(Lx/nx)*i;
pE[i+j*ny].nd[2].y=(Ly/ny)*(j+1);
pE[i+j*ny].nd[2].number=i+(nx+1)*(j+1); //NO.2
pE[i+j*ny].nd[3].x=(Lx/nx)*(i+1);
pE[i+j*ny].nd[3].y=(Ly/ny)*(j+1);
pE[i+j*ny].nd[3].number=i+1+(nx+1)*(j+1); //NO.3
}
//System.out.println("计算基函数系数值...\n");
int j=1,k=2;
for(int id=0;id<nx*ny;id++)
{
pE[id].A=(pE[id].nd[1].x-pE[id].nd[0].x)*(pE[id].nd[2].y-pE[id].nd[0].y);
D=pE[id].A;
}
//--------------------------------------计算单元有限元特征式系数矩阵可以不再分两个三角形循环
// for(int i=0;i<nx*ny;i++)
// {
// pE[i].strainRatioX=((pE[i].nd[0].velocityX-pE[i].nd[2].velocityX)*(pE[i].nd[1].y-pE[i].nd[3].y)-(pE[i].nd[1].velocityX-pE[i].nd[3].velocityX)*(pE[i].nd[0].y-pE[i].nd[2].y))/(2*pE[i].A);
// pE[i].strainRatioY=-((pE[i].nd[0].velocityY-pE[i].nd[2].velocityY)*(pE[i].nd[1].x-pE[i].nd[3].x)-(pE[i].nd[1].velocityY-pE[i].nd[3].velocityY)*(pE[i].nd[0].x-pE[i].nd[2].x))/(2*pE[i].A);
// pE[i].strainRatioXY=((pE[i].nd[0].velocityY-pE[i].nd[2].velocityY)*(pE[i].nd[1].y-pE[i].nd[3].y)-(pE[i].nd[1].velocityY-pE[i].nd[3].velocityY+(pE[i].nd[0].velocityX-pE[i].nd[2].velocityX)*(pE[i].nd[1].x-pE[i].nd[3].x)-(pE[i].nd[1].velocityX-pE[i].nd[3].velocityX)*(pE[i].nd[0].x-pE[i].nd[2].x))*(pE[i].nd[0].y-pE[i].nd[2].y))/(2*pE[i].A);
//
// pE[i].strainX=((pE[i].nd[0].displacementX-pE[i].nd[2].displacementX)*(pE[i].nd[1].y-pE[i].nd[3].y)-(pE[i].nd[1].displacementX-pE[i].nd[3].displacementX)*(pE[i].nd[0].y-pE[i].nd[2].y))/(2*pE[i].A);;
// pE[i].strainY=-((pE[i].nd[0].displacementY-pE[i].nd[2].displacementY)*(pE[i].nd[1].x-pE[i].nd[3].x)-(pE[i].nd[1].displacementY-pE[i].nd[3].displacementY)*(pE[i].nd[0].x-pE[i].nd[2].x))/(2*pE[i].A);;
//
//
// pE[i].strainXY=((pE[i].nd[0].displacementY-pE[i].nd[2].displacementY)*(pE[i].nd[1].y-pE[i].nd[3].y)-(pE[i].nd[1].displacementY-pE[i].nd[3].displacementY+(pE[i].nd[0].displacementX-pE[i].nd[2].displacementX)*(pE[i].nd[1].x-pE[i].nd[3].x)-(pE[i].nd[1].displacementX-pE[i].nd[3].displacementX)*(pE[i].nd[0].x-pE[i].nd[2].x))*(pE[i].nd[0].y-pE[i].nd[2].y))/(2*pE[i].A);
//
// //平面应变
// //pE[i].p=E*(pE[i].strainX+pE[i].strainY)*(1-u)/(3*(1-2*u));
// //pE[i].Sx=E*((2-u)*pE[i].strainX-(1+u)*pE[i].strainY)/(3*(1+u));
// //pE[i].Sy=E*((2-u)*pE[i].strainY-(1+u)*pE[i].strainX)/(3*(1+u));
//
// //平面应力
// pE[i].p=E*(pE[i].strainX+pE[i].strainY)/(3*(1-2*u));
// pE[i].Sx=E*(2*pE[i].strainX-pE[i].strainY)/(3*(1+u));
// pE[i].Sy=E*(2*pE[i].strainY-pE[i].strainX)/(3*(1+u));
// pE[i].sigmaX=pE[i].p+pE[i].Sx;
// pE[i].sigmaY=pE[i].p+pE[i].Sy;
// pE[i].sigmaXY=E*pE[i].strainXY/(1+u);
//
// //节点x方向力
// FDMElementTriangle_2D tmpA=null,tmpB=null,tmpC=null,tmpD=null,tmpE=null,tmpF=null,tmpG=null,tmpH=null,tmpI=null,tmpJ=null;
// if((i-nx-1)<0)
// tmpA=(FDMElementTriangle_2D)pE[i].clone();
// else{
// tmpA=pE[i-nx-1];
// }
//
// if((i-nx)<0)
// tmpB=(FDMElementTriangle_2D)pE[i].clone();
// else{
// tmpB=pE[i-nx];
// }
//
// if((i-1)<0)
// tmpC=(FDMElementTriangle_2D)pE[i].clone();
// else{
// tmpC=pE[i-1];
// }
// Fx[pE[i].nd[0].number ]+=(tmpA.sigmaX*(tmpA.nd[1].y-tmpA.nd[2].y)-tmpA.sigmaXY*(tmpA.nd[1].x-tmpA.nd[2].x)+
// tmpB.sigmaX*(tmpB.nd[3].y-tmpB.nd[0].y)-tmpB.sigmaXY*(tmpB.nd[3].x-tmpB.nd[0].x)+
// pE[i].sigmaX*(pE[i].nd[2].y-pE[i].nd[1].y)-pE[i].sigmaXY*(pE[i].nd[2].x-pE[i].nd[1].x)+
// tmpC.sigmaX*(tmpC.nd[0].y-tmpC.nd[3].y)-tmpC.sigmaXY*(tmpC.nd[0].x-tmpC.nd[3].x)
// );
//
// if((i-nx+1)<0){
// tmpE=(FDMElementTriangle_2D)pE[i].clone();
// }
// else{
// tmpE=pE[i-nx+1];
// }
// if((i+1)>=nx*ny){
// tmpD=(FDMElementTriangle_2D)pE[i].clone();
// }
// else{
// tmpD=pE[i+1];
// }
// Fx[pE[i].nd[1].number ]+=(tmpB.sigmaX*(tmpB.nd[1].y-tmpB.nd[2].y)-tmpB.sigmaXY*(tmpB.nd[1].x-tmpB.nd[2].x)+
// tmpE.sigmaX*(tmpE.nd[3].y-tmpE.nd[0].y)-tmpE.sigmaXY*(tmpE.nd[3].x-tmpE.nd[0].x)+
// tmpD.sigmaX*(tmpD.nd[2].y-tmpD.nd[1].y)-tmpD.sigmaXY*(tmpD.nd[2].x-tmpD.nd[1].x)+
// pE[i].sigmaX*(pE[i].nd[0].y-pE[i].nd[3].y)-pE[i].sigmaXY*(pE[i].nd[0].x-pE[i].nd[3].x)
// );
// if((i-1)<0){
// tmpF=(FDMElementTriangle_2D)pE[i].clone();
// }
// else{
// tmpF=pE[i-1];
// }
// if((i+nx)>=nx*ny){
// tmpG=(FDMElementTriangle_2D)pE[i].clone();
// }
// else{
// tmpG=pE[i+nx];
// }
// if((i+1)>=nx*ny){
// tmpH=(FDMElementTriangle_2D)pE[i].clone();
// }
// else{
// tmpH=pE[i+1];
// }
// if((i+nx-1)>=nx*ny){
// tmpI=(FDMElementTriangle_2D)pE[i].clone();
// }
// else{
// tmpI=pE[i+nx-1];
// }
// if((i+nx+1)>=nx*ny){
// tmpJ=(FDMElementTriangle_2D)pE[i].clone();
// }
// else{
// tmpJ=pE[i+nx+1];
// }
// Fx[pE[i].nd[2].number ]+=(tmpF.sigmaX*(tmpF.nd[1].y-tmpF.nd[2].y)-tmpF.sigmaXY*(tmpF.nd[1].x-tmpF.nd[2].x)+
// pE[i].sigmaX*(pE[i].nd[3].y-pE[i].nd[0].y)-pE[i].sigmaXY*(pE[i].nd[3].x-pE[i].nd[0].x)+
// tmpG.sigmaX*(tmpG.nd[2].y-tmpG.nd[1].y)-tmpG.sigmaXY*(tmpG.nd[2].x-tmpG.nd[1].x)+
// tmpI.sigmaX*(tmpI.nd[0].y-tmpI.nd[3].y)-tmpI.sigmaXY*(tmpI.nd[0].x-tmpI.nd[3].x)
// );
// Fx[pE[i].nd[3].number ]+=(pE[i].sigmaX*(pE[i].nd[1].y-pE[i].nd[2].y)-pE[i].sigmaXY*(pE[i].nd[1].x-pE[i].nd[2].x)+
// tmpH.sigmaX*(tmpH.nd[3].y-tmpH.nd[0].y)-tmpH.sigmaXY*(tmpH.nd[3].x-tmpH.nd[0].x)+
// tmpJ.sigmaX*(tmpJ.nd[2].y-tmpJ.nd[1].y)-tmpJ.sigmaXY*(tmpJ.nd[2].x-tmpJ.nd[1].x)+
// tmpG.sigmaX*(tmpG.nd[0].y-tmpG.nd[3].y)-tmpG.sigmaXY*(tmpG.nd[0].x-tmpG.nd[3].x)
// );
//
// //节点y方向力
// FDMElementTriangle_2D tmpyA=null,tmpyB=null,tmpyC=null,tmpyE=null,tmpyF=null,tmpyG=null,tmpyH=null,tmpyI=null,tmpyJ=null;
// if((i-nx-1)<0)
// tmpyA=(FDMElementTriangle_2D)pE[i].clone();
// else{
// tmpyA=pE[i-nx-1];
// }
//
// if((i-nx)<0)
// tmpyB=(FDMElementTriangle_2D)pE[i].clone();
// else{
// tmpyB=pE[i-nx];
// }
//
// if((i-1)<0)
// tmpyC=(FDMElementTriangle_2D)pE[i].clone();
// else{
// tmpyC=pE[i-1];
// }
// Fy[pE[i].nd[0].number ]+=(tmpyA.sigmaXY*(tmpyA.nd[1].y-tmpyA.nd[2].y)-tmpyA.sigmaY*(tmpyA.nd[1].x-tmpyA.nd[2].x)+
// tmpyB.sigmaXY*(tmpyB.nd[3].y-tmpyB.nd[0].y)-tmpyB.sigmaY*(tmpyB.nd[3].x-tmpyB.nd[0].x)+
// pE[i].sigmaXY*(pE[i].nd[2].y-pE[i].nd[1].y)-pE[i].sigmaY*(pE[i].nd[2].x-pE[i].nd[1].x)+
// tmpyC.sigmaXY*(tmpyC.nd[0].y-tmpyC.nd[3].y)-tmpyC.sigmaY*(tmpyC.nd[0].x-tmpyC.nd[3].x)
// );
//
// if((i-nx+1)<0){
// tmpyE=(FDMElementTriangle_2D)pE[i].clone();
// }
// else{
// tmpyE=pE[i-nx+1];
// }
//
// //
// if((i+1)>=nx*ny){
// tmpyG=(FDMElementTriangle_2D)pE[i].clone();
// }
// else{
// tmpyG=pE[i+1];
// }
// if((i+nx)>=nx*ny){
// tmpyH=(FDMElementTriangle_2D)pE[i].clone();
// }
// else{
// tmpyH=pE[i+nx];
// }
// if((i+nx-1)>=nx*ny){
// tmpyI=(FDMElementTriangle_2D)pE[i].clone();
// }
// else{
// tmpyI=pE[i+nx-1];
// }
// if((i+nx+1)>=nx*ny){
// tmpyJ=(FDMElementTriangle_2D)pE[i].clone();
// }
// else{
// tmpyJ=pE[i+nx+1];
// }
// Fy[pE[i].nd[1].number ]+=(tmpyB.sigmaXY*(tmpyB.nd[1].y-tmpyB.nd[2].y)-tmpyB.sigmaY*(tmpyB.nd[1].x-tmpyB.nd[2].x)+
// tmpyE.sigmaXY*(tmpyE.nd[3].y-tmpyE.nd[0].y)-tmpyE.sigmaY*(tmpyE.nd[3].x-tmpyE.nd[0].x)+
// tmpyG.sigmaXY*(tmpyG.nd[2].y-tmpyG.nd[1].y)-tmpyG.sigmaY*(tmpyG.nd[2].x-tmpyG.nd[1].x)+
// pE[i].sigmaXY*(pE[i].nd[0].y-pE[i].nd[3].y)-pE[i].sigmaY*(pE[i].nd[0].x-pE[i].nd[3].x)
// );
// if((i-1)<0){
// tmpyF=(FDMElementTriangle_2D)pE[i].clone();
// }
// else{
// tmpyF=pE[i-1];
// }
// Fy[pE[i].nd[2].number ]+=(tmpyF.sigmaXY*(tmpyF.nd[1].y-tmpyF.nd[2].y)-tmpyF.sigmaY*(tmpyF.nd[1].x-tmpyF.nd[2].x)+
// pE[i].sigmaXY*(pE[i].nd[3].y-pE[i].nd[0].y)-pE[i].sigmaY*(pE[i].nd[3].x-pE[i].nd[0].x)+
// tmpyH.sigmaXY*(tmpyH.nd[2].y-tmpyH.nd[1].y)-tmpyH.sigmaY*(tmpyH.nd[2].x-tmpyH.nd[1].x)+
// tmpyI.sigmaXY*(tmpyI.nd[0].y-tmpyI.nd[3].y)-tmpyI.sigmaY*(tmpyI.nd[0].x-tmpyI.nd[3].x)
// );
// Fy[pE[i].nd[3].number ]+=(pE[i].sigmaXY*(pE[i].nd[1].y-pE[i].nd[2].y)-pE[i].sigmaY*(pE[i].nd[1].x-pE[i].nd[2].x)+
// tmpyG.sigmaXY*(tmpyG.nd[3].y-tmpyG.nd[0].y)-tmpyG.sigmaY*(tmpyG.nd[3].x-tmpyG.nd[0].x)+
// tmpyJ.sigmaXY*(tmpyJ.nd[2].y-tmpyJ.nd[1].y)-tmpyJ.sigmaY*(tmpyJ.nd[2].x-tmpyJ.nd[1].x)+
// tmpyH.sigmaXY*(tmpyH.nd[0].y-tmpyH.nd[3].y)-tmpyH.sigmaY*(tmpyH.nd[0].x-tmpyH.nd[3].x)
// );
//
//
// total_accX[pE[i].nd[0].number]=Fx[pE[i].nd[0].number ]/(rou*(Lx/nx));
// total_accX[pE[i].nd[0].number]=Fx[pE[i].nd[0].number ]/(rou*(Ly/ny));
//
// total_accX[pE[i].nd[1].number]=Fx[pE[i].nd[1].number ]/(rou*(Lx/nx));
// total_accX[pE[i].nd[1].number]=Fx[pE[i].nd[1].number ]/(rou*(Ly/ny));
//
// total_accX[pE[i].nd[2].number]=Fx[pE[i].nd[2].number ]/(rou*(Lx/nx));
// total_accX[pE[i].nd[2].number]=Fx[pE[i].nd[2].number ]/(rou*(Ly/ny));
//
// total_accX[pE[i].nd[3].number]=Fx[pE[i].nd[3].number ]/(rou*(Lx/nx));
// total_accX[pE[i].nd[3].number]=Fx[pE[i].nd[3].number ]/(rou*(Ly/ny));
// }
//---------------------------计算三角形单元f函数向量不用多重积分公式,直接代用单元三角形面积的三分之一
int idx=0;
for(int i=0;i<nx*ny;i++)
{ for(idx=0;idx<3;idx++)
{
pE[i].fi[idx]=(Lx/nx)*(Ly/ny); // Respaired
}
}
//单元矩阵元素累加到总体矩阵相应的位置上,共四个节点，两个单元
// 2 |-------|3
// | I . |
// | . II |
// 0 |-------|1
//double dBig=pow(10,20); //边界条件对角线扩大法处理所用的大数
// double Ur=0.0; //边界条件1(边界条件2通过Calerkin弱解表达式自动满足)
double Ur=1.23 ; //外力，牛
/*
*
* jimmy 上下拉伸约束
*/
//上端2,3节点，Y方向拉力 索引为5,7
// for(int i=0;i<nx+1;i++)
// {
//
// Fx[(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(int n=0;n<nx+1;n++)
{
total_accX[n]=0;
total_accY[n]=0;
total_velocityX[n]=0;
total_velocityY[n]=0;
total_displacementX[n]=0;
total_displacementY[n]=0;
}
//第一步求t=0加速度的，结果在Fx中
//假设t=0，速度和位移为0
// for(int i=0;i<2*iNode;i++)
// tmp[i]=Fx[i];
//假设初始加速度为确定值
for(int i=iNode-(nx+1);i<iNode;i++)
{
total_displacementY[i]=0.01;
total_accY[i]=-1;
Fy[i]=Math.sqrt(E)*total_velocityY[i];
}
//增加时间步长的内容
// double p=0,Sx=0,Sy=0,Sxy=0;
double p=0,Sx=0,Sy=0,Sxy=0;
for(int step=1;step<T/dataT;step++)
{
for(int n=0;n<nx+1;n++)
{
total_accX[n]=0;
total_accY[n]=0;
total_velocityX[n]=0;
total_velocityY[n]=0;
total_displacementX[n]=0;
total_displacementY[n]=0;
}
//double tmpX[]=total_velocityX;
// total_strainRatioX[i]=(-total_velocityX[i])/(Lx);
//FDMElementTriangle_2D tmpA=null,tmpB=null,tmpC=null,tmpE=null,tmpF=null;
//节点y方向力
//FDMElementTriangle_2D tmpyA=null,tmpyB=null,tmpyC=null,tmpyE=null,tmpyF=null;
for(int i=0;i<nx*ny;i++)
{
pE[i].nd[0].velocityX=total_velocityX[pE[i].nd[0].number];
pE[i].nd[0].velocityY=total_velocityY[pE[i].nd[0].number];
pE[i].nd[1].velocityX=total_velocityX[pE[i].nd[1].number];
pE[i].nd[1].velocityY=total_velocityY[pE[i].nd[1].number];
pE[i].nd[2].velocityX=total_velocityX[pE[i].nd[2].number];
pE[i].nd[2].velocityY=total_velocityY[pE[i].nd[2].number];
pE[i].nd[3].velocityX=total_velocityX[pE[i].nd[3].number];
pE[i].nd[3].velocityY=total_velocityY[pE[i].nd[3].number];
pE[i].nd[0].displacementX=total_displacementX[pE[i].nd[0].number];
pE[i].nd[0].displacementY=total_displacementY[pE[i].nd[0].number];
pE[i].nd[1].displacementX=total_displacementX[pE[i].nd[1].number];
pE[i].nd[1].displacementY=total_displacementY[pE[i].nd[1].number];
pE[i].nd[2].displacementX=total_displacementX[pE[i].nd[2].number];
pE[i].nd[2].displacementY=total_displacementY[pE[i].nd[2].number];
pE[i].nd[3].displacementX=total_displacementX[pE[i].nd[3].number];
pE[i].nd[3].displacementY=total_displacementY[pE[i].nd[3].number];
pE[i].nd[0].accX=total_accX[pE[i].nd[0].number];
pE[i].nd[0].accY=total_accY[pE[i].nd[0].number];
pE[i].nd[1].accX=total_accX[pE[i].nd[1].number];
pE[i].nd[1].accY=total_accY[pE[i].nd[1].number];
pE[i].nd[2].accX=total_accX[pE[i].nd[2].number];
pE[i].nd[2].accY=total_accY[pE[i].nd[2].number];
pE[i].nd[3].accX=total_accX[pE[i].nd[3].number];
pE[i].nd[3].accY=total_accY[pE[i].nd[3].number];
}
for(int i=0;i<nx*ny;i++)
{
pE[i].strainRatioX=((pE[i].nd[0].velocityX-pE[i].nd[2].velocityX)*(pE[i].nd[1].y-pE[i].nd[3].y)-(pE[i].nd[1].velocityX-pE[i].nd[3].velocityX)*(pE[i].nd[0].y-pE[i].nd[2].y))/(2*pE[i].A);
pE[i].strainRatioY=-((pE[i].nd[0].velocityY-pE[i].nd[2].velocityY)*(pE[i].nd[1].x-pE[i].nd[3].x)-(pE[i].nd[1].velocityY-pE[i].nd[3].velocityY)*(pE[i].nd[0].x-pE[i].nd[2].x))/(2*pE[i].A);
pE[i].strainRatioXY=((pE[i].nd[0].velocityY-pE[i].nd[2].velocityY)*(pE[i].nd[1].y-pE[i].nd[3].y)-(pE[i].nd[1].velocityY-pE[i].nd[3].velocityY+(pE[i].nd[0].velocityX-pE[i].nd[2].velocityX)*(pE[i].nd[1].x-pE[i].nd[3].x)-(pE[i].nd[1].velocityX-pE[i].nd[3].velocityX)*(pE[i].nd[0].x-pE[i].nd[2].x))*(pE[i].nd[0].y-pE[i].nd[2].y))/(2*pE[i].A);
pE[i].strainX=((pE[i].nd[0].displacementX-pE[i].nd[2].displacementX)*(pE[i].nd[1].y-pE[i].nd[3].y)-(pE[i].nd[1].displacementX-pE[i].nd[3].displacementX)*(pE[i].nd[0].y-pE[i].nd[2].y))/(2*pE[i].A);;
pE[i].strainY=-((pE[i].nd[0].displacementY-pE[i].nd[2].displacementY)*(pE[i].nd[1].x-pE[i].nd[3].x)-(pE[i].nd[1].displacementY-pE[i].nd[3].displacementY)*(pE[i].nd[0].x-pE[i].nd[2].x))/(2*pE[i].A);;
pE[i].strainXY=((pE[i].nd[0].displacementY-pE[i].nd[2].displacementY)*(pE[i].nd[1].y-pE[i].nd[3].y)-(pE[i].nd[1].displacementY-pE[i].nd[3].displacementY+(pE[i].nd[0].displacementX-pE[i].nd[2].displacementX)*(pE[i].nd[1].x-pE[i].nd[3].x)-(pE[i].nd[1].displacementX-pE[i].nd[3].displacementX)*(pE[i].nd[0].x-pE[i].nd[2].x))*(pE[i].nd[0].y-pE[i].nd[2].y))/(2*pE[i].A);
//平面应变strainZ=0
//pE[i].p=E*(pE[i].strainX+pE[i].strainY)/(3*(1-2*u));
//pE[i].Sx+=E*((2-3*u)*pE[i].strainRatioX/(1-2*u)-(1-3*u)*pE[i].strainY/(1-2*u))/(3*(1+u));
//pE[i].Sy+=E*((2-3*u)*pE[i].strainRatioY/(1-2*u)-(1-3*u)*pE[i].strainX/(1-2*u)))/(3*(1+u));
//pE[i].Sz+=-E*(1+2*u)(pE[i].strainRatioX+pE[i].strainRatioX)/(3*(1-u+2*u*u));
pE[i].p=pE[i].p+dataT*E*(pE[i].strainRatioX+pE[i].strainRatioX)/(3*(1-2*u));
pE[i].Sx=pE[i].Sx+dataT*E*(2*pE[i].strainRatioX-pE[i].strainRatioY)/(3*(1+u))+0.5*(pE[i].Sx-pE[i].Sy)*(Math.cos(2*dataT*pE[i].strainRatioXY)-1)-pE[i].Sxy*Math.sin(2*dataT*pE[i].strainRatioXY);
pE[i].Sy=pE[i].Sy+dataT*E*(2*pE[i].strainRatioY-pE[i].strainRatioX)/(3*(1+u))-0.5*(pE[i].Sx-pE[i].Sy)*(Math.cos(2*dataT*pE[i].strainRatioXY)-1)+pE[i].Sxy*Math.sin(2*dataT*pE[i].strainRatioXY);;
pE[i].Sxy=pE[i].Sxy+dataT*E*(pE[i].strainRatioXY)/((1+u))+0.5*(pE[i].Sx-pE[i].Sy)*(Math.sin(2*dataT*pE[i].strainRatioXY))+pE[i].Sxy*(Math.cos(2*dataT*pE[i].strainRatioXY)-1);
pE[i].sigmaX=pE[i].p+pE[i].Sx;
pE[i].sigmaY=pE[i].p+pE[i].Sy;
pE[i].sigmaXY=E*pE[i].Sxy;
// pE[i].p=E*(pE[i].strainX+pE[i].strainY)/(3*(1-2*u));
// pE[i].Sx=E*((2-3*u)*pE[i].strainRatioX/(1-2*u)-(1-3*u)*pE[i].strainY/(1-2*u))/(3*(1+u));
// pE[i].Sy=E*((2-3*u)*pE[i].strainRatioY/(1-2*u)-(1-3*u)*pE[i].strainX/(1-2*u))/(3*(1+u));
//
//
// pE[i].sigmaX=pE[i].p+pE[i].Sx;
// pE[i].sigmaY=pE[i].p+pE[i].Sy;
//pE[i].sigmaXY=E*pE[i].Sxy;
//平面应力
//只有y方向的位移，应变，（以及引起的x的方向的）
//p=(E*displacementX/Lx/3);
//Sy=(E*2*strainRatioY*dataT/3/2);
//Sy=(E/(1+u)*(strainRatioY-0.33333*total_eratio)*dataT);//-0.5*(Sx-Sy)*(Math.cos(2*dataT*strainRatioXY)-1)-Sxy*Math.sin(2*dataT*strainRatioXY)
// Sxy=(E/(1+u)*(strainRatioXY)*dataT);
//
// pE[i].p=E*(pE[i].strainX+pE[i].strainY)*(1-u)/(3*(1-2*u));
// pE[i].Sx=pE[i].Sx+dataT*E*((2-u)*pE[i].strainRatioX-(1+u)*pE[i].strainRatioY)/(3*(1+u))+0.5*(pE[i].Sx-pE[i].Sy)*(Math.cos(2*dataT*pE[i].strainRatioXY)-1)-pE[i].Sxy*Math.sin(2*dataT*pE[i].strainRatioXY);
//
// pE[i].Sy=pE[i].Sy+dataT*E*((2-u)*pE[i].strainRatioY-(1+u)*pE[i].strainRatioX)/(3*(1+u))-0.5*(pE[i].Sx-pE[i].Sy)*(Math.cos(2*dataT*pE[i].strainRatioXY)-1)+pE[i].Sxy*Math.sin(2*dataT*pE[i].strainRatioXY);;
//
// pE[i].Sxy=pE[i].Sxy+dataT*E*(pE[i].strainRatioXY)/((1+u))+0.5*(pE[i].Sx-pE[i].Sy)*(Math.sin(2*dataT*pE[i].strainRatioXY))+pE[i].Sxy*(Math.cos(2*dataT*pE[i].strainRatioXY)-1);
//
// pE[i].sigmaX=pE[i].p+pE[i].Sx;
// pE[i].sigmaY=pE[i].p+pE[i].Sy;
// pE[i].sigmaXY=E*pE[i].Sxy;
//节点x方向力
FDMElementTriangle_2D tmpA=null,tmpB=null,tmpC=null,tmpD=null,tmpE=null,tmpF=null,tmpG=null,tmpH=null,tmpI=null,tmpJ=null;
if((i-nx-1)<0)
tmpA=(FDMElementTriangle_2D)pE[i].clone();
else{
tmpA=pE[i-nx-1];
}
if((i-nx)<0)
tmpB=(FDMElementTriangle_2D)pE[i].clone();
else{
tmpB=pE[i-nx];
}
if((i-1)<0)
tmpC=(FDMElementTriangle_2D)pE[i].clone();
else{
tmpC=pE[i-1];
}
Fx[pE[i].nd[0].number ]=(tmpA.sigmaX*(tmpA.nd[1].y-tmpA.nd[2].y)-tmpA.sigmaXY*(tmpA.nd[1].x-tmpA.nd[2].x)+
tmpB.sigmaX*(tmpB.nd[3].y-tmpB.nd[0].y)-tmpB.sigmaXY*(tmpB.nd[3].x-tmpB.nd[0].x)+
pE[i].sigmaX*(pE[i].nd[2].y-pE[i].nd[1].y)-pE[i].sigmaXY*(pE[i].nd[2].x-pE[i].nd[1].x)+
tmpC.sigmaX*(tmpC.nd[0].y-tmpC.nd[3].y)-tmpC.sigmaXY*(tmpC.nd[0].x-tmpC.nd[3].x)
);
if((i-nx+1)<0){
tmpE=(FDMElementTriangle_2D)pE[i].clone();
}
else{
tmpE=pE[i-nx+1];
}
if((i+1)>=nx*ny){
tmpD=(FDMElementTriangle_2D)pE[i].clone();
}
else{
tmpD=pE[i+1];
}
Fx[pE[i].nd[1].number ]=(tmpB.sigmaX*(tmpB.nd[1].y-tmpB.nd[2].y)-tmpB.sigmaXY*(tmpB.nd[1].x-tmpB.nd[2].x)+
tmpE.sigmaX*(tmpE.nd[3].y-tmpE.nd[0].y)-tmpE.sigmaXY*(tmpE.nd[3].x-tmpE.nd[0].x)+
tmpD.sigmaX*(tmpD.nd[2].y-tmpD.nd[1].y)-tmpD.sigmaXY*(tmpD.nd[2].x-tmpD.nd[1].x)+
pE[i].sigmaX*(pE[i].nd[0].y-pE[i].nd[3].y)-pE[i].sigmaXY*(pE[i].nd[0].x-pE[i].nd[3].x)
);
if((i-1)<0){
tmpF=(FDMElementTriangle_2D)pE[i].clone();
}
else{
tmpF=pE[i-1];
}
if((i+nx)>=nx*ny){
tmpG=(FDMElementTriangle_2D)pE[i].clone();
}
else{
tmpG=pE[i+nx];
}
if((i+1)>=nx*ny){
tmpH=(FDMElementTriangle_2D)pE[i].clone();
}
else{
tmpH=pE[i+1];
}
if((i+nx-1)>=nx*ny){
tmpI=(FDMElementTriangle_2D)pE[i].clone();
}
else{
tmpI=pE[i+nx-1];
}
if((i+nx+1)>=nx*ny){
tmpJ=(FDMElementTriangle_2D)pE[i].clone();
}
else{
tmpJ=pE[i+nx+1];
}
Fx[pE[i].nd[2].number ]=(tmpF.sigmaX*(tmpF.nd[1].y-tmpF.nd[2].y)-tmpF.sigmaXY*(tmpF.nd[1].x-tmpF.nd[2].x)+
pE[i].sigmaX*(pE[i].nd[3].y-pE[i].nd[0].y)-pE[i].sigmaXY*(pE[i].nd[3].x-pE[i].nd[0].x)+
tmpG.sigmaX*(tmpG.nd[2].y-tmpG.nd[1].y)-tmpG.sigmaXY*(tmpG.nd[2].x-tmpG.nd[1].x)+
tmpI.sigmaX*(tmpI.nd[0].y-tmpI.nd[3].y)-tmpI.sigmaXY*(tmpI.nd[0].x-tmpI.nd[3].x)
);
Fx[pE[i].nd[3].number ]=(pE[i].sigmaX*(pE[i].nd[1].y-pE[i].nd[2].y)-pE[i].sigmaXY*(pE[i].nd[1].x-pE[i].nd[2].x)+
tmpH.sigmaX*(tmpH.nd[3].y-tmpH.nd[0].y)-tmpH.sigmaXY*(tmpH.nd[3].x-tmpH.nd[0].x)+
tmpJ.sigmaX*(tmpJ.nd[2].y-tmpJ.nd[1].y)-tmpJ.sigmaXY*(tmpJ.nd[2].x-tmpJ.nd[1].x)+
tmpG.sigmaX*(tmpG.nd[0].y-tmpG.nd[3].y)-tmpG.sigmaXY*(tmpG.nd[0].x-tmpG.nd[3].x)
);
//节点y方向力
FDMElementTriangle_2D tmpyA=null,tmpyB=null,tmpyC=null,tmpyE=null,tmpyF=null,tmpyG=null,tmpyH=null,tmpyI=null,tmpyJ=null;
if((i-nx-1)<0)
tmpyA=(FDMElementTriangle_2D)pE[i].clone();
//tmpyA=pE[i-nx-1].
else{
tmpyA=pE[i-nx-1];
}
if((i-nx)<0)
tmpyB=(FDMElementTriangle_2D)pE[i].clone();
else{
tmpyB=pE[i-nx];
}
if((i-1)<0)
tmpyC=(FDMElementTriangle_2D)pE[i].clone();
else{
tmpyC=pE[i-1];
}
Fy[pE[i].nd[0].number ]=(tmpyA.sigmaXY*(tmpyA.nd[1].y-tmpyA.nd[2].y)-tmpyA.sigmaY*(tmpyA.nd[1].x-tmpyA.nd[2].x)+
tmpyB.sigmaXY*(tmpyB.nd[3].y-tmpyB.nd[0].y)-tmpyB.sigmaY*(tmpyB.nd[3].x-tmpyB.nd[0].x)+
pE[i].sigmaXY*(pE[i].nd[2].y-pE[i].nd[1].y)-pE[i].sigmaY*(pE[i].nd[2].x-pE[i].nd[1].x)+
tmpyC.sigmaXY*(tmpyC.nd[0].y-tmpyC.nd[3].y)-tmpyC.sigmaY*(tmpyC.nd[0].x-tmpyC.nd[3].x)
);
if((i-nx+1)<0){
tmpyE=(FDMElementTriangle_2D)pE[i].clone();
}
else{
tmpyE=pE[i-nx+1];
}
//
if((i+1)>=nx*ny){
tmpyG=(FDMElementTriangle_2D)pE[i].clone();
}
else{
tmpyG=pE[i+1];
}
if((i+nx)>=nx*ny){
tmpyH=(FDMElementTriangle_2D)pE[i].clone();
}
else{
tmpyH=pE[i+nx];
}
if((i+nx-1)>=nx*ny){
tmpyI=(FDMElementTriangle_2D)pE[i].clone();
}
else{
tmpyI=pE[i+nx-1];
}
if((i+nx+1)>=nx*ny){
tmpyJ=(FDMElementTriangle_2D)pE[i].clone();
}
else{
tmpyJ=pE[i+nx+1];
}
Fy[pE[i].nd[1].number ]=(tmpyB.sigmaXY*(tmpyB.nd[1].y-tmpyB.nd[2].y)-tmpyB.sigmaY*(tmpyB.nd[1].x-tmpyB.nd[2].x)+
tmpyE.sigmaXY*(tmpyE.nd[3].y-tmpyE.nd[0].y)-tmpyE.sigmaY*(tmpyE.nd[3].x-tmpyE.nd[0].x)+
tmpyG.sigmaXY*(tmpyG.nd[2].y-tmpyG.nd[1].y)-tmpyG.sigmaY*(tmpyG.nd[2].x-tmpyG.nd[1].x)+
pE[i].sigmaXY*(pE[i].nd[0].y-pE[i].nd[3].y)-pE[i].sigmaY*(pE[i].nd[0].x-pE[i].nd[3].x)
);
if((i-1)<0){
tmpyF=(FDMElementTriangle_2D)pE[i].clone();
}
else{
tmpyF=pE[i-1];
}
Fy[pE[i].nd[2].number ]=(tmpyF.sigmaXY*(tmpyF.nd[1].y-tmpyF.nd[2].y)-tmpyF.sigmaY*(tmpyF.nd[1].x-tmpyF.nd[2].x)+
pE[i].sigmaXY*(pE[i].nd[3].y-pE[i].nd[0].y)-pE[i].sigmaY*(pE[i].nd[3].x-pE[i].nd[0].x)+
tmpyH.sigmaXY*(tmpyH.nd[2].y-tmpyH.nd[1].y)-tmpyH.sigmaY*(tmpyH.nd[2].x-tmpyH.nd[1].x)+
tmpyI.sigmaXY*(tmpyI.nd[0].y-tmpyI.nd[3].y)-tmpyI.sigmaY*(tmpyI.nd[0].x-tmpyI.nd[3].x)
);
Fy[pE[i].nd[3].number ]=(pE[i].sigmaXY*(pE[i].nd[1].y-pE[i].nd[2].y)-pE[i].sigmaY*(pE[i].nd[1].x-pE[i].nd[2].x)+
tmpyG.sigmaXY*(tmpyG.nd[3].y-tmpyG.nd[0].y)-tmpyG.sigmaY*(tmpyG.nd[3].x-tmpyG.nd[0].x)+
tmpyJ.sigmaXY*(tmpyJ.nd[2].y-tmpyJ.nd[1].y)-tmpyJ.sigmaY*(tmpyJ.nd[2].x-tmpyJ.nd[1].x)+
tmpyH.sigmaXY*(tmpyH.nd[0].y-tmpyH.nd[3].y)-tmpyH.sigmaY*(tmpyH.nd[0].x-tmpyH.nd[3].x)
);
// total_accX[pE[i].nd[0].number]+=Fx[pE[i].nd[0].number ]/(rou*(Lx/nx));
// total_accY[pE[i].nd[0].number]+=Fy[pE[i].nd[0].number ]/(rou*(Ly/ny));
//
// total_accX[pE[i].nd[1].number]+=Fx[pE[i].nd[1].number ]/(rou*(Lx/nx));
// total_accY[pE[i].nd[1].number]+=Fy[pE[i].nd[1].number ]/(rou*(Ly/ny));
//
// total_accX[pE[i].nd[2].number]+=Fx[pE[i].nd[2].number ]/(rou*(Lx/nx));
// total_accY[pE[i].nd[2].number]+=Fy[pE[i].nd[2].number ]/(rou*(Ly/ny));
//
// total_accX[pE[i].nd[3].number]+=Fx[pE[i].nd[3].number ]/(rou*(Lx/nx));
// total_accY[pE[i].nd[3].number]+=Fy[pE[i].nd[3].number ]/(rou*(Ly/ny));
tmpA=null;
tmpB=null;
tmpC=null;
tmpE=null;
tmpF=null;
tmpyA=null;
tmpyB=null;
tmpyC=null;
tmpyE=null;
tmpyF=null;
}
for(int ii=0;ii<iNode;ii++)
{
total_accX[ii]=Fx[ii]/(rou*(Lx*Ly/(nx*ny)))/2;
total_accY[ii]=Fy[ii]/(rou*(Lx*Ly/(ny*ny)))/2;
total_velocityX[ii]=total_velocityX[ii]+total_accX[ii]*dataT;
total_velocityY[ii]=total_velocityY[ii]+total_accY[ii]*dataT;
total_displacementX[ii]=total_displacementX[ii]+total_velocityX[ii]*dataT;
total_displacementY[ii]=total_displacementY[ii]+total_velocityY[ii]*dataT;
//if(step%3==0)
System.out.println("节点编号："+3+",y方向位移 :"+displacementX+"\t, y方向位移应力: "+(p)+"y方向加速度： "+acc);
//Gauss(total_accX,Fx,iNode);
//total_accX[i]=-(p+Sy)/(rou*Lx);
}
if(step%8==0)
System.out.println(total_displacementY[iNode-1]);
System.out.println("节点编号："+iNode+"\t"+Fx[2*iNode-2]+"\t"+Fx[2*iNode-1]);
System.out.println(Fx[2*iNode-2]+"\t"+Fx[2*iNode-1]);
}
}
catch(Exception e)
{
e.printStackTrace();
System.out.println("Error occured...\n");
}
return ;
}
}
1842
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
