|
|
【答案】应助回帖
#include "stdafx.h"
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#define PI 3.14159265
#define MaxGrid 750
double p[MaxGrid+5][MaxGrid+5],T[MaxGrid+5][MaxGrid+5],c[MaxGrid+5][MaxGrid+5],r[MaxGrid+5][MaxGrid+5],difT,difTx,difTy;
double hp[MaxGrid+5][MaxGrid+5],hpp[MaxGrid+5][MaxGrid+5];
//=====================初试化相场和温度场==========================
void IninPhaseAndTemp()
{
int i,j;
for(i=0;i<MaxGrid+2;i++)
{
for(j=0;j<MaxGrid+2;j++)
{
if(i+j<26)//256
{
p[j]=0;
T[j]=1574;// 1574
c[j]=0.3994;//;0.026842
}
else
{
p[j]=1;
T[j]=1574;//////// 1574
c[j]=0.40831;//0.07214
}
}
}
for(i=0;i<MaxGrid+4;i++)
{
p[1]=p[3];
p[MaxGrid+2]=p[MaxGrid];
p[1]=p[3];
p[MaxGrid+2]=p[MaxGrid];
T[1]=T[3];
T[MaxGrid+2]=T[MaxGrid];
T[1]=T[3];
T[MaxGrid+2]=T[MaxGrid];
c[0]=c[4];
c[1]=c[3];
c[MaxGrid+3]=c[MaxGrid-1];
c[MaxGrid+2]=c[MaxGrid];
c[0]=c[4];
c[1]=c[3];
c[MaxGrid+3]=c[MaxGrid-1];
c[MaxGrid+2]=c[MaxGrid];
}
}
//=====================求解相场==========================
void Phase()
{
double difpx,difpy,difpxx,difpxy,difpyy,angx,ang,angy,difcx,difcy,difcxx,difcyy,difcxxx,difcyyy,difcxxy,difcyyx,difcxxxx,difcyyyy,difcxxyy,difcyyxx;
double Pa,Pb,Ca,Cb,Cc,Cd,Gp,hppp,Gpp,Gppp,Ha,Hb,Haa,Hbb,Dc,Dcc;//,Cb
///////计算相场和温度场///////
for(j=2;j<MaxGrid+2;j++)
{
if(p[2][j]<1.0 && p[2][j+1]==1)
{
MinGrid=j+1;
printf("\n MinGrid=%d p[2][j]=%f p[2][j+1]=%f\n",j,p[2][j],p[2][j+1]);
break;
}
}
/////////////////////////////Noise////////////////////////////
for(i=2;i<MinGrid+2;i++)
{
for(j=2;j<MinGrid+2;j++)
{
r[j]=rand()/32767.0*2.0-1;
}
}
//////////////////////////////////////////////////////////////////
for(i=2;i<MinGrid+2;i++)
{
for(j=2;j<MinGrid+2;j++)
{
////////////////////////////////////////////////////////////////
difpx=(p[i+1][j]-p[i-1][j])/(2*deltx);
difpy=(p[j+1]-p[j-1])/(2*delty);
difpxy=(p[i+1][j+1]-p[i-1][j+1]-p[i+1][j-1]+p[i-1][j-1])/(4*deltx*delty);
difpxx=(p[i+1][j]-2*p[j]+p[i-1][j])/(deltx*deltx);
difpyy=(p[j+1]-2*p[j]+p[j-1])/(delty*delty);
difcx=(c[i+1][j]-c[i-1][j])/(2*deltx);
difcy=(c[j+1]-c[j-1])/(2*delty);
difcxx=(c[i+1][j]-2*c[j]+c[i-1][j])/(deltx*deltx);
difcyy=(c[j+1]-2*c[j]+c[j-1])/(delty*delty);
difcxxx=(c[i+2][j]-2*c[i+1][j]+2*c[i-1][j]-c[i-2][j])/(2*deltx*deltx*deltx);
difcyyy=(c[j+2]-2*c[j+1]+2*c[j-1]-c[j-2])/(2*delty*delty*delty);
difcxxy=(c[i+1][j+1]-2*c[j+1]+c[i-1][j+1]-c[i+1][j-1]+2*c[j-1]-c[i-1][j-1])/(2*deltx*deltx*deltx);
difcyyx=(c[i+1][j+1]-2*c[i+1][j]+c[i+1][j-1]-c[i-1][j+1]+2*c[i-1][j]-c[i-1][j-1])/(2*delty*delty*delty);
difcxxyy=(c[i+1][j+1]-2*c[j+1]+c[i-1][j+1]+c[i+1][j-1]-2*c[j-1]+c[i-1][j-1]-2*c[i+1][j]-2*c[i-1][j]+4*c[j])/(deltx*deltx*deltx*deltx);
difcyyxx=(c[i+1][j+1]-2*c[i+1][j]+c[i+1][j-1]+c[i-1][j+1]-2*c[i-1][j]+c[i-1][j-1]-2*c[j+1]-2*c[j-1]+4*c[j])/(deltx*deltx*deltx*deltx);
difcxxxx=(c[i+2][j]-4*c[i+1][j]+6*c[j]-4*c[i-1][j]+c[i-2][j])/(deltx*deltx*deltx*deltx);
difcyyyy=(c[j+2]-4*c[j+1]+6*c[j]-4*c[j-1]+c[j-2])/(delty*delty*delty*delty);
difTx=(T[i+1][j]-T[i-1][j])/(2*deltx);
difTy=(T[j+1]-T[j-1])/(2*delty);
difT=(T[i+1][j]+T[i-1][j]+T[j-1]+T[j+1]-4.0*T[j])/(deltx*deltx);
/////////////////////////////////////////////////////////////////
Gp=(p[j]*p[j])*(1-p[j])*(1-p[j]);
Gpp=(4*p[j]*p[j]*p[j]-6*p[j]*p[j]+2*p[j]);
Gppp=(12*p[j]*p[j]-12*p[j]+2);
hp[j]=(10*p[j]*p[j]*p[j]-15*p[j]*p[j]*p[j]*p[j]+6*p[j]*p[j]*p[j]*p[j]*p[j]);
hpp[j]=(30*p[j]*p[j]*p[j]*p[j]-60*p[j]*p[j]*p[j]+30*p[j]*p[j]);
hppp=(120*p[j]*p[j]*p[j]-180*p[j]*p[j]+60*p[j]);
Ha=Wa*Gpp+hpp[j]*La*(1/T[j]-1/Tma);
Hb=Wb*Gpp+hpp[j]*Lb*(1/T[j]-1/Tmb);
Haa=Wa*Gppp+hppp*La*(1/T[j]-1/Tma);
Hbb=Wb*Gppp+hppp*Lb*(1/T[j]-1/Tmb);
Dc=Ds+hp[j]*(Dl-Ds);
Dcc=hpp[j]*(Dl-Ds);
//////////////////////////////////////////////////////////////////
if(difpx!=0)
{
ang=atan(difpy/difpx);
}
else if(difpy==0)
{
ang=PI/4;
}
else
{
ang=PI/2;
}
if(difpx==0 && difpy==0)
{
angx=1;
angy=1;
}
else
{
angx=(difpx*difpxy-difpy*difpxx)/(difpx*difpx+difpy*difpy);
angy=(difpx*difpyy-difpy*difpxy)/(difpx*difpx+difpy*difpy);
}
////////////////p/////////////////////
Pa=ea*ea*(16*gamma*(cos(4*ang)+gamma*cos(8*ang))*(angx*difpy-angy*difpx)
-(angx*difpx+angy*difpy)*(8*gamma*sin(4*ang)+4*gamma*gamma*sin(8*ang))
+(difpxx+difpyy)*(1+2*gamma*cos(4*ang)+gamma*gamma*cos(4*ang)*cos(4*ang)));//
Pb=9*r[j]*16*Gp*((1-c[j])*Ha+c[j]*Hb);
Ca=Vm/R*((1-2*c[j])*Dc*(Hb-Ha)*(difcx*difpx+difcy*difpy)
+c[j]*(1-c[j])*Dcc*(Hb-Ha)*(difpx*difpx+difpy*difpy)
+c[j]*(1-c[j])*Dc*(Hbb-Haa)*(difpx*difpx+difpy*difpy)
+c[j]*(1-c[j])*Dc*(Hb-Ha)*(difpxx+difpyy)
);
Cb=Dcc*(difpx*difcx+difpy*difcy)+Dc*(difcxx+difcyy);
Cc=2*Vm/R*(La-Lb)*((1-2*c[j])*Dc*hp[j]*1/T[j]*1/T[j]*(difcx*difTx+difcy*difTy)
+c[j]*(1-c[j])*Dcc*hp[j]*1/T[j]*1/T[j]*(difpx*difTx+difpy*difTy)
+c[j]*(1-c[j])*Dc*hpp[j]*1/T[j]*1/T[j]*(difpx*difTx+difpy*difTy)
-2*c[j]*(1-c[j])*Dc*hp[j]*1/T[j]*1/T[j]*1/T[j]*(difTx*difTx+difTy*difTy)
+c[j]*(1-c[j])*Dc*hp[j]*1/T[j]*1/T[j]*difT);
Cd=-Vm/R*del*del*((1-2*c[j])*Dc*(difcx*difcxxx+difcx*difcyyx+difcy*difcxxy+difcy*difcyyy)
+c[j]*(1-c[j])*Dcc*(difpx*difcxxx+difpx*difcyyx+difpy*difcxxy+difpy*difcyyy)
+c[j]*(1-c[j])*Dc*(difcxxxx+difcyyxx+difcxxyy+difcyyyy));
newp[j]=p[j]+delt*((1-c[j])*Ma+c[j]*Mb)*(Pa-(1-c[j])*Ha-c[j]*Hb-Pb);
newc[j]=c[j]+delt*(Ca+Cb+Cc+Cd);
}
}
//////////////////////////////////////////////////////
if((num%2==0))
{
for(i=2;i<MinGrid+2;i++)
{
for(j=2;j<MinGrid+2;j++)
{
DD[j]=T[j-1]+BB*T[j]+T[j+1]-(deltx*deltx)*((1-c[j])*La+c[j]*Lb)*hpp[j]*(newp[j]-p[j])/(DT*Cp*delt)
-(deltx*deltx)*(La-Lb)*(1-hp[j])*(newc[j]-c[j])/(DT*Cp*delt);
}
}
for(j=2;j<MinGrid+2;j++)
{
Cm[MinGrid+1][j]=1.0/(B-1.0);
Cn[MinGrid+1][j]=DD[MinGrid+1][j]/(B-1.0);
for(i=MinGrid;i>1;i--)//求系数Cm,Cn;
{
Cm[j]=1.0/(B-Cm[i+1][j]);
Cn[j]=(Cn[i+1][j]+DD[j])/(B-Cm[i+1][j]);
}
newT[2][j]=(Cn[2][j]+DD[2][j])/(-Cm[2][j]+B-1.0);
for(i=3;i<MinGrid+2;i++)
{
newT[j]=Cm[j]*newT[i-1][j]+Cn[j];
}
}
}
///////////////////////////////////////y direction//////////////////////////////////////////////////////////
else
{
for(i=2;i<MinGrid+2;i++)
{
for(j=2;j<MinGrid+2;j++)
{
DD[j]=T[i-1][j]+BB*T[j]+T[i+1][j]-(deltx*deltx)*((1-c[j])*La+c[j]*Lb)*hpp[j]*(newp[j]-p[j])/(DT*Cp*delt)
-(deltx*deltx)*(La-Lb)*(1-hp[j])*(newc[j]-c[j])/(DT*Cp*delt);
}
}
for(i=2;i<MinGrid+2;i++)
{
Cm[MinGrid+1]=1.0/(B-1.0);//有问题
Cn[MinGrid+1]=DD[MinGrid+1]/(B-1.0);//
for(j=MinGrid;j>1;j--)//求系数Cm,Cn;
{
Cm[j]=1.0/(B-Cm[j+1]);
Cn[j]=(Cn[j+1]+DD[j])/(B-Cm[j+1]);
}
newT[2]=(Cn[2]+DD[2])/(-Cm[2]+B-1.0);
for(j=3;j<MinGrid+2;j++)
{
newT[j]=Cm[j]*newT[j-1]+Cn[j];
}
}
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////
for(i=2;i<MinGrid+2;i++) //MinGrid
{
for(j=2;j<MinGrid+2;j++) //MinGrid
{
p[j]=newp[j];
c[j]=newc[j];
// T[j]=newT[j];
}
}
///////////////设置边界条件//////////
for(i=0;i<MaxGrid+4;i++)
{
p[0]=p[4];
p[1]=p[3];
p[MaxGrid+3]=p[MaxGrid-1];
p[MaxGrid+2]=p[MaxGrid];
p[0]=p[4];
p[1]=p[3];
p[MaxGrid+3]=p[MaxGrid-1];
p[MaxGrid+2]=p[MaxGrid];
T[1]=T[3]; | |
回复此楼
