进行结构响应分析的Wilson法
进行结构响应分析的Wilson法，是直接积分的基本算法之一，相信进行有限元动力分析编程的诸位都能用上（这是vc代码）void CGlobalElement::WilsonMethod()
{
int loop,loop1,iBuf;
int nTimeStep,nLoadTimeStep,nGroundAccData,nTotalDOF,nFreeDOF;
double *adGroundAcc;
double dBuf,dBuf1,dAlfa,dPsi,da,da1,da2,da3,da4,da5;
CSparseMatrix smatMass,smatDK;
CMatrix matDisp,matVec,matAcc,matBuf,matBuf1;
ofstream fout;
ifstream fin;
if(m_apEle.GetSize()==0) return;
m_bFlagDynamicAnalyzed=true;
nTotalDOF=m_Node.GetTotalDOF();
nFreeDOF=m_Node.GetFreeDOF();
matDisp.Realloc(nFreeDOF,1);
matVec.Realloc(nFreeDOF,1);
matAcc.Realloc(nFreeDOF,1);
matBuf.Realloc(nFreeDOF,1);
matBuf1.Realloc(nFreeDOF,1);
GetMassMatrix(smatMass);
m_smatGK.ElementBufRealloc();
m_smatGK=0.0;
for(loop=0;loop<m_nEle;loop++){
m_apEle->StiffAssemble(m_smatGK);
}
dBuf=(m_adEigen*m_adEigen-m_adEigen*m_adEigen);
dAlfa=2.0*m_adEigen*m_adEigen*(m_adDampRatio*m_adEigen
-m_adDampRatio*m_adEigen)/dBuf;
dPsi=2.0*(m_adDampRatio*m_adEigen-m_adDampRatio*m_adEigen)/dBuf;
dBuf1=m_dWilsonXita*m_dTimeStep;
dBuf=1.0/dBuf1;
da=3.0*dBuf*(2.0*dBuf+dAlfa);
da1=1.0+3.0*dBuf*dPsi;
da2=3.0+dBuf1*dAlfa/2.0;
da3=dBuf1*dPsi/2.0;
da4=6.0*dBuf+3.0*dAlfa;
da5=3.0*dPsi;
nTimeStep=(int)(m_dResponseDuration/m_dTimeStep);
nLoadTimeStep=(int)(m_dLoadDuration/m_dTimeStep);
smatDK=m_smatGK;
smatDK*=da1;
smatMass*=da;
smatDK+=smatMass;
dBuf=1/da;
smatMass*=dBuf;
matDisp=0.0;matVec=0.0;matAcc=0.0;
fout.open("dydata.tmp",ios::binary);
if(m_iDynamicLoadType==DYNAMIC_INSTANT_LOAD){
m_iCurLoadGroup=0;
GetLoadVector();
for(loop=0;loop<nTimeStep;loop++){
matBuf=matVec*da4+matAcc*da2;
matBuf=smatMass*matBuf;
matBuf1=matVec*da5+matAcc*da3;
matBuf1=m_smatGK*matBuf1;
for(loop1=0;loop1<nFreeDOF;loop1++)
m_adDisp=matBuf(loop1,0)+matBuf1(loop1,0);
if(loop==0){
for(loop1=0;loop1<nFreeDOF;loop1++)
m_adDisp+=m_dWilsonXita*m_adLoadVector;
}
else if(loop==nLoadTimeStep){
for(loop1=0;loop1<nFreeDOF;loop1++)
m_adDisp-=m_dWilsonXita*m_adLoadVector;
}
if(!smatDK.LdltSolve(nFreeDOF,m_adDisp)) return;
for(loop1=0;loop1<nFreeDOF;loop1++)
matBuf(loop1,0)=m_adDisp;
dBuf=6.0/(m_dWilsonXita*m_dWilsonXita*m_dWilsonXita*m_dTimeStep*m_dTimeStep);
dBuf1=-6.0/(m_dWilsonXita*m_dWilsonXita*m_dTimeStep);
matBuf=matBuf*dBuf+matVec*dBuf1+matAcc*(1.0-3.0/m_dWilsonXita);
matDisp+=matVec*m_dTimeStep+(matBuf+matAcc*2.0)*(m_dTimeStep*m_dTimeStep/6.0);
matVec+=(matBuf+matAcc)*(m_dTimeStep/2.0);
matAcc=matBuf;
for(loop1=0;loop1<nFreeDOF;loop1++)
m_adDisp=matDisp(loop1,0);
fout.write((char*)m_adDisp,sizeof(double)*nFreeDOF);
}
}
else{
fin.open(m_sGroundAccFile,ios::nocreate);
if(!fin.good()){
AfxMessageBox("Can't open seismic acceleration file", MB_OK, 0 );
return;
}
dBuf1=0.0;
nGroundAccData=0;
while(!fin.eof()){
fin>>dBuf;
nGroundAccData++;
if(dBuf1<fabs(dBuf)) dBuf1=fabs(dBuf);
}
fin.close();
dBuf=m_dPeakAcc/dBuf1;
adGroundAcc=new double ;
fin.open(m_sGroundAccFile,ios::nocreate);
adGroundAcc=0.0;
for(loop=1;loop<nGroundAccData;loop++){
fin>>adGroundAcc;
adGroundAcc*=dBuf;
}
fin.close();
for(loop=0;loop<nTimeStep;loop++){
matBuf=matVec*da4+matAcc*da2;
if(loop<nGroundAccData){
for(loop1=0;loop1<m_nNode;loop1++){
iBuf=m_Node.GetXDOFIndex(loop1);
if(iBuf<nFreeDOF){
matBuf(iBuf,0)+=m_dWilsonXita*(adGroundAcc-adGroundAcc);
}
}
}
matBuf=smatMass*matBuf;
matBuf1=matVec*da5+matAcc*da3;
matBuf1=m_smatGK*matBuf1;
for(loop1=0;loop1<nFreeDOF;loop1++)
m_adDisp=matBuf(loop1,0)+matBuf1(loop1,0);
if(!smatDK.LdltSolve(nFreeDOF,m_adDisp)) return;
for(loop1=0;loop1<nFreeDOF;loop1++)
matBuf(loop1,0)=m_adDisp;
dBuf=6.0/(m_dWilsonXita*m_dWilsonXita*m_dWilsonXita*m_dTimeStep*m_dTimeStep);
dBuf1=-6.0/(m_dWilsonXita*m_dWilsonXita*m_dTimeStep);
matBuf=matBuf*dBuf+matVec*dBuf1+matAcc*(1.0-3.0/m_dWilsonXita);
matDisp+=matVec*m_dTimeStep+(matBuf+matAcc*2.0)*(m_dTimeStep*m_dTimeStep/6.0);
matVec+=(matBuf+matAcc)*(m_dTimeStep/2.0);
matAcc=matBuf;
for(loop1=0;loop1<nFreeDOF;loop1++)
m_adDisp=matDisp(loop1,0);
fout.write((char*)m_adDisp,sizeof(double)*nFreeDOF);
}
delete adGroundAcc;
}
fout.close();
} 正需要啊，
谢谢
页:
[1]