24小时热门版块排行榜    

CyRhmU.jpeg
查看: 1893  |  回复: 4

浩淼323

金虫 (小有名气)

[求助] 求助pudn账号程序下载已有2人参与

回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐,对您同样有帮助:

已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

浩淼323

金虫 (小有名气)

各位大神,帮个忙,谢谢啦!!!!
2楼2014-01-03 22:00:26
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

fish.yfyh

铜虫 (小有名气)

【答案】应助回帖

感谢参与,应助指数 +1
3楼2014-01-03 22:10:51
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

qobitroway

铁虫 (初入文坛)

【答案】应助回帖

#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];   
        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 output_tempresult()   
{   
    int i,j;  char filename[30];      
    FILE *fp;   
   
    if(num>50&&num%5000==0)   
    {   
        
        sprintf(filename,"%d.plt",num);   
        fp=fopen(filename,"w";   
      
        fprintf(fp,"%s","title='pu.txt'\n";   
        //fprintf(fp,"%s","variables= i , j , p , c\n";   
        fprintf(fp,"%s","variables= i , j , p , t, c\n";   
        fprintf(fp,"%s","zone I=750,J=750,F=point\n";   
    //  for(i=2;i<MaxGrid+2;i++)   
         for(j=2;j<MaxGrid+2;j++)   
        {   
    ///     for(j=2;j<MaxGrid+2;j++)   
            for(i=2;i<MaxGrid+2;i++)   
            {   
              fprintf(fp,"%d  %d ",i,j);//   
              fprintf(fp,"%3.5f   %4.4f   %7.5f\n",p[j],T[j],c[j]);   
              //fprintf(fp,"%3.2f  %7.5f\n", p[j], c[j]);   
              //fprintf(fp,"%3.2f   %3.2f\n",p[j],T[j]);   
            }   
            fprintf(fp,"\n";              
        }   
            fprintf(fp,"\n";      
             fclose(fp);   
      }   
    }   
   
int main(int argc, char* argv[])   
{   
            
    t=0;   
    num=0;     
    int numv=0;   
    int kkp=10,kkr,nold=0,n=0;   
    double vv;//,pxr,pyr,rv,   
    FILE *fpv;   
    char filename1[30];   
    sprintf(filename1,"%d.txt",111);   
//////////////////////////////////////////   
   printf("The programming is running,not to shut down!";   
   
   printf("f\n",Ma);   
   //printf("i=%f\n   j=%f\n   p=%15.14f\n",Mf,W,eps);   
   IninPhaseAndTemp();      
/////////////////////////////////////////      
   while(t<2000)   
   {      
    Phase();   
  //for(j=2;j<MaxGrid+2;j++)   
  //{   
    for(i=2;i<MaxGrid+2;i++)   
    {            
       if(p[0]*p[i+1][0]<=0.8)   
       {            
           kkr=i;   
           if((kkr-kkp)==5)   
           {   
               for(i=0;i<MinGrid+2;i++)     //MinGrid   
               {   
                 for(j=0;j<MinGrid+2;j++) //MinGrid   
                 {     
           
                 if(T[j]>Txy)  Txy=T[j] ;   
                 if(p[j]<0.9)  gu=gu+1 ;   
   
                 }   
               }     
             vv=5*deltx/(((n-nold)*delt));   
              
             fpv=fopen(filename1,"a+";            
   
             fprintf(fpv,"%d ",kkr);   
             fprintf(fpv,"%f  %f  %f\n",vv,Txy,gu/562500);   
             fclose(fpv);   
             kkp=kkr;   
             Txy=1574;gu=0;   
             nold=n;   
           }   
       }   
    }     
    output_tempresult();   
    t=t+delt;   
    n=n+1;     
    num++;     
    printf("num=%d",num);   
   }     
    return 0;   
}
4楼2014-01-10 09:31:11
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

hetuotuo

新虫 (初入文坛)

引用回帖:
3楼: Originally posted by fish.yfyh at 2014-01-03 22:10:51
源代码可以免费得到啊
http://read.pudn.com/downloads188/sourcecode/database/885766/Phase.cpp__.htm

求告知,如何免费获得,地址是怎么得到的
5楼2016-01-01 22:31:47
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 浩淼323 的主题更新
信息提示
请填处理意见