www.pudn.com > mainstructure.rar > mainstructure.f90, change:2015-01-04,size:15962b


 
    module global 
     
    implicit none 
     
    integer                  :: nnodes,ncells,nedges,nwing 
    !for node number,cell number,edge number,node number of wing respectively 
     
    integer,allocatable      :: iedge(:,:),icell(:,:),celledge(:),wing(:,:) 
    !iedge(1,i)      :   edge's start point 'a' 
    !iedge(2,i)      :   edge's end point 'b' 
    !iedge(3,i)      :   edge's left cell  (all positive) 
    !iedge(4,i)      :   positive for right cell, -1 for wall boundary, -2 for farfield boundary 
    !icell           :   triangle cell's three nodes index 
    !celledge(ncells):   the number of a cell's edges which are not boundaries 
    !wing(2,nwing)   :   for wing's uper and lower nodes respectively 
     
    real*8,allocatable       :: xy(:,:),dxy(:,:),vol(:) 
    !xy(2,nedges)    :   cartesian coordinates xy(1,i) for x, xy(2,i) for y 
    !dxy(2,nedges)   :   discrepence of edge's two point, dxy(1,i) for x, dxy(2,i) for y 
    !vol(ncells)     :   cell's volume(area in 2d) 
     
    real*8,parameter         ::gama=1.4,rouinf=1.223,pinf=101325.,k2=0.8,k4=0.02 
    real*8::machinf,uinf,vinf,cinf,angle,cfl,epsilon 
    common machinf,uinf,vinf,cinf,angle,cfl,epsilon 
    !for mach number,,P,u,v,local speed of sound of the freeflow respectively 
    !angle for angle of attack 
    !epsilon for smoothing coefficient 
     
    real*8,allocatable       ::w(:,:),q(:,:),d(:,:),scale_a(:),delta_t(:),r(:,:) 
    real*8                   ::r_restrain(4) 
    !scale_a(ncells) :   for scaling factor and t 
    !delta_t(ncells) :   time step 
    !r(4,ncells)     :   residual error 
    !r_restrain(4)   :   restrain criterion 
     
    end module 
     
    program main 
     
    use global 
    implicit none 
     
    call readgrid 
    print*,'Reading grid data complete!' 
     
    call initialize 
    print*,'Reading initial data successful!' 
     
    call iterate 
    call results 
     
    pause 
     
    end program 
     
    subroutine  readgrid 
     
    use global 
    implicit none 
     
    integer::i,j,k,celltype 
     
    open(10,file='naca0012-structure.grd') 
    read(10,*)  celltype,nnodes,nedges,ncells 
     
    allocate(iedge(4,nedges)) 
    allocate(xy(2,nnodes)) 
    allocate(dxy(2,nedges)) 
    allocate(icell(4,ncells)) 
    allocate(celledge(ncells)) 
    allocate(vol(ncells)) 
    allocate(w(4,ncells)) 
    allocate(q(4,ncells)) 
    allocate(d(4,ncells)) 
    allocate(delta_t(ncells)) 
    allocate(scale_a(ncells)) 
    allocate(r(4,ncells)) 
     
    nwing=0 
    celledge=0 
     
    do i=1,nnodes 
        read(10,*)   xy(1,i),xy(2,i) 
    end do 
     
    do i=1,nedges 
        read(10,*) iedge(1,i),iedge(2,i),iedge(3,i),iedge(4,i) 
    end do 
     
    do i=1,ncells 
        read(10,*) icell(1,i),icell(2,i),icell(3,i),icell(4,i) 
    end do 
     
    do i=1,ncells 
        read(10,*) vol(i) 
    end do 
     
    close(10) 
     
    do i=1,nedges 
         
        dxy(1,i)=xy(1,iedge(2,i))-xy(1,iedge(1,i)) 
        dxy(2,i)=xy(2,iedge(2,i))-xy(2,iedge(1,i)) 
         
        if(iedge(4,i)>0) then 
             
            celledge(iedge(3,i))=celledge(iedge(3,i))+1 
            celledge(iedge(4,i))=celledge(iedge(4,i))+1 
             
        end if 
         
        if(iedge(4,i)==-1) nwing=nwing+1 
         
    end do 
     
    nwing=nwing/2+1 
    allocate(wing(2,nwing)) 
    j=0 
    k=0 
     
    do i=1,nedges 
         
        if(iedge(4,i)==-1) then 
             
            if(xy(2,iedge(1,i))>=0.) then 
                 
                j=j+1 
                wing(1,j)=iedge(1,i) 
                 
            else if(xy(2,iedge(1,i))<0.) then 
                 
                k=k+1 
                wing(2,k)=iedge(1,i) 
                 
            end if 
             
        end if 
         
    end do 
     
    do i=1,nwing-1 
        do j=i+1,nwing 
             
            if(xy(1,wing(1,j))<xy(1,wing(1,i))) then 
                 
                k=wing(1,i) 
                wing(1,i)=wing(1,j) 
                wing(1,j)=k 
                 
            end if 
             
            if(j<nwing-1) then 
                if(xy(1,wing(2,j))<xy(1,wing(2,i))) then 
                     
                    k=wing(2,i) 
                    wing(2,i)=wing(2,j) 
                    wing(2,j)=k 
                     
                end if 
            end if 
             
        end do 
    end do 
     
    do i=nwing-1,2,-1 
         
        wing(2,i)=wing(2,i-1) 
         
    end do 
     
    wing(2,1)=wing(1,1) 
    wing(2,nwing)=wing(1,nwing) 
     
    return 
     
    end subroutine 
     
    subroutine initialize 
     
    use global 
    implicit none 
     
    real*8::pi=4.*atan(1.) 
     
    print*,'Please check the existence of the file "initial.dat" whose content are:' 
    print*,'Mach    Angle of attack(degree)    cfl    epsilon' 
    pause 
     
    open(10,file='initial.dat') 
    read(10,*) machinf,angle,cfl,epsilon 
    close(10) 
     
    angle=angle*pi/180. 
    cinf=sqrt(pinf*gama/rouinf) 
    uinf=machinf*cinf*cos(angle) 
    vinf=machinf*cinf*sin(angle) 
     
    w(1,:)=rouinf 
    w(2,:)=rouinf*uinf 
    w(3,:)=rouinf*vinf 
    w(4,:)=pinf/(gama-1.)+0.5*rouinf*(uinf**2+vinf**2) 
     
    return 
     
    end subroutine 
     
    subroutine iterate 
     
    use global 
    implicit none 
     
    integer::i 
    real*8::r0(4) 
     
    i=0 
    r_restrain=1.0 
     
    open(20,file='restrain.dat') 
    write(*,"(2X,A5,2X,A7,6X,A1,11X,A1,8X,A7/)") 'Times','Density','U','V','Energy' 
     
    do while(log10(r_restrain(1))>-6) 
         
        i=i+1 
        call runge_kutta 
        if(i==1) r0=r_restrain 
        r_restrain=r_restrain/r0 
         
        if(mod(i,100)==0) write(*,"(I6,2X,F7.3,3X,F7.3,6X,F7.3,6X,F7.3)") i,log10(r_restrain(:)) 
         
        write(20,"(I6,2X,F7.3,3X,F7.3,6X,F7.3,6X,F7.3)") i,log10(r_restrain(:)) 
         
    end do 
     
    close(20) 
    return 
     
    end subroutine 
     
    subroutine runge_kutta 
     
    use global 
    implicit none 
     
    integer::i,j 
    real*8::w0(4,ncells),alpha(4) 
     
    w0=w 
    alpha(1)=1.0/4.0 
    alpha(2)=1.0/3.0 
    alpha(3)=1.0/2.0 
    alpha(4)=1.0 
    r_restrain=0. 
     
    do i=1,4 
         
        call compute_flux 
         
        if(i==1) then 
             
            call artificial_d 
             
        end if 
         
        call r_smooth 
         
        do j=1,ncells 
             
            w(:,j)=w0(:,j)+alpha(i)*delta_t(j)*r(:,j) 
             
        end do 
         
    end do 
     
    do i=1,ncells 
         
        r_restrain(:)=r_restrain(:)+r(:,i)**2 
         
    end do 
     
    r_restrain=sqrt(r_restrain/ncells) 
     
    end subroutine 
     
    subroutine compute_flux 
     
    use global 
    implicit none 
     
    integer::i,j,k,p 
    real*8::axy,wi(4),pi,roui,ui,vi,ci,zi,flux(4) 
    real*8::u0,v0,c0,rou0,p0,mach0,vninf,vtinf,vn0,vt0,vn,vt,c 
     
    q=0. 
    scale_a=0. 
     
    do i=1,nedges 
         
        k=iedge(3,i) 
        p=iedge(4,i) 
        axy=sqrt(dxy(1,i)**2+dxy(2,i)**2) 
         
        if(p>0) then 
             
            wi(:)=0.5*(w(:,k)+w(:,p)) 
            roui=wi(1) 
            ui=wi(2)/wi(1) 
            vi=wi(3)/wi(1) 
            pi=(gama-1.)*(wi(4)-0.5*roui*(ui**2+vi**2)) 
            ci=sqrt(pi*gama/roui) 
            zi=ui*dxy(2,i)-vi*dxy(1,i) 
            flux(1)=zi*wi(1) 
            flux(2)=zi*wi(2)+pi*dxy(2,i) 
            flux(3)=zi*wi(3)-pi*dxy(1,i) 
            flux(4)=zi*(wi(4)+pi) 
             
        else if(p==-1) then 
             
            wi(:)=w(:,k) 
            roui=wi(1) 
            ui=wi(2)/wi(1) 
            vi=wi(3)/wi(1) 
            pi=(gama-1.)*(wi(4)-0.5*roui*(ui**2+vi**2)) 
            ci=sqrt(pi*gama/roui) 
            zi=0. 
            flux(1)=0. 
            flux(2)=pi*dxy(2,i) 
            flux(3)=-pi*dxy(1,i) 
            flux(4)=0. 
             
        else if(p==-2) then 
             
            rou0=w(1,k) 
            u0=w(2,k)/w(1,k) 
            v0=w(3,k)/w(1,k) 
            p0=(gama-1.)*(w(4,k)-0.5*rou0*(u0**2+v0**2)) 
            c0=sqrt(p0*gama/rou0) 
            vninf=(uinf*dxy(2,i)-vinf*dxy(1,i))/axy 
            vtinf=(uinf*dxy(1,i)+vinf*dxy(2,i))/axy 
            vn0=(u0*dxy(2,i)-v0*dxy(1,i))/axy 
            vt0=(u0*dxy(1,i)+v0*dxy(2,i))/axy 
            vn=0.5*(vninf+vn0) 
            c=0.5*(cinf+c0) 
            mach0=abs(vn/c) 
             
            if(mach0<=1.0) then 
                if(vn<=0.0) then     !subsonic in-flow 
                     
                     
                    vt=vtinf 
                    ui=(vn*dxy(2,i)+vt*dxy(1,i))/axy 
                    vi=(vt*dxy(2,i)-vn*dxy(1,i))/axy 
                    ci=c 
                    roui=rouinf*((ci**2)/(cinf**2))**(1./(gama-1.)) 
                    pi=roui*(ci**2)/gama 
                    wi(1)=roui 
                    wi(2)=roui*ui 
                    wi(3)=roui*vi 
                    wi(4)=pi/(gama-1.)+0.5*roui*(ui**2+vi**2) 
                    zi=vn*axy 
                    flux(1)=zi*wi(1) 
                    flux(2)=zi*wi(2)+pi*dxy(2,i) 
                    flux(3)=zi*wi(3)-pi*dxy(1,i) 
                    flux(4)=zi*(wi(4)+pi) 
                     
                else                  !subsonic out-flow 
                     
                    vt=vt0 
                    ui=(vn*dxy(2,i)+vt*dxy(1,i))/axy 
                    vi=(vt*dxy(2,i)-vn*dxy(1,i))/axy 
                    ci=c 
                    roui=rou0*((ci**2)/(c0**2))**(1./(gama-1.)) 
                    pi=roui*ci**2/gama 
                    wi(1)=roui 
                    wi(2)=roui*ui 
                    wi(3)=roui*vi 
                    wi(4)=pi/(gama-1.)+0.5*roui*(ui**2+vi**2) 
                    zi=vn*axy 
                    flux(1)=zi*wi(1) 
                    flux(2)=zi*wi(2)+pi*dxy(2,i) 
                    flux(3)=zi*wi(3)-pi*dxy(1,i) 
                    flux(4)=zi*(wi(4)+pi) 
                 
                end if 
                 
            else 
                vn=vninf 
                if(vn<=0.0) then     !supersonic in-flow 
                     
                    ui=uinf 
                    vi=vinf 
                    ci=cinf 
                    pi=pinf 
                    roui=rouinf 
                    wi(1)=roui 
                    wi(2)=roui*ui 
                    wi(3)=roui*vi 
                    wi(4)=pi/(gama-1.)+0.5*roui*(ui**2+vi**2) 
                    zi=vn*axy 
                    flux(1)=zi*wi(1) 
                    flux(2)=zi*wi(2)+pi*dxy(2,i) 
                    flux(3)=zi*wi(3)-pi*dxy(1,i) 
                    flux(4)=zi*(wi(4)+pi) 
                     
                else                  !supersonic out-flow 
                     
                    ui=u0 
                    vi=v0 
                    ci=c0 
                    pi=p0 
                    roui=rou0 
                    wi(1)=roui 
                    wi(2)=roui*ui 
                    wi(3)=roui*vi 
                    wi(4)=pi/(gama-1.)+0.5*roui*(ui**2+vi**2) 
                    zi=ui*dxy(2,i)-vi*dxy(1,i) 
                    flux(1)=zi*wi(1) 
                    flux(2)=zi*wi(2)+pi*dxy(2,i) 
                    flux(3)=zi*wi(3)-pi*dxy(1,i) 
                    flux(4)=zi*(wi(4)+pi) 
                     
                end if 
                 
            end if 
             
        else 
             
            print*,'Data error!' 
            pause 
            stop 
             
        end if 
         
        q(:,k)=q(:,k)+flux(:) 
        scale_a(k)=scale_a(k)+abs(zi)+ci*axy 
         
        if(p>0) then 
            q(:,p)=q(:,p)-flux(:) 
            scale_a(p)=scale_a(p)+abs(zi)+ci*axy 
        end if 
         
    end do 
     
    return 
     
    end subroutine 
     
    subroutine artificial_d 
     
    use global 
    implicit none 
     
    integer::i,k,p 
    real*8::shock_v,epson2,epson4,w_laplace(4,ncells),d2(4),d4(4) 
    real*8::u,v,c,pp,pk 
     
    w_laplace=0. 
    d=0. 
     
    do i=1,nedges 
         
        k=iedge(3,i) 
        p=iedge(4,i) 
         
        if(p>0) then 
             
            w_laplace(:,k)=w_laplace(:,k)+0.5*(w(:,p)-w(:,k)) 
            w_laplace(:,p)=w_laplace(:,p)+0.5*(w(:,k)-w(:,p)) 
             
        end if 
         
    end do 
     
    do i=1,nedges 
         
        k=iedge(3,i) 
        p=iedge(4,i) 
         
        if(p>0) then 
             
            pp=(gama-1.)*(w(4,p)-0.5*(w(2,p)**2+w(3,p)**2)/w(1,p)) 
            pk=(gama-1.)*(w(4,k)-0.5*(w(2,k)**2+w(3,k)**2)/w(1,k)) 
            shock_v=abs((pp-pk)/(pp+pk)) 
            epson2=k2*shock_v 
            epson4=max(0.0,(k4-epson2)) 
            d2(:)=0.5*(scale_a(k)+scale_a(p))*epson2*(w(:,p)-w(:,k)) 
            d4(:)=0.5*(scale_a(k)+scale_a(p))*epson4*(w_laplace(:,k)-w_laplace(:,p)) 
            d(:,k)=d(:,k)+d2(:)+d4(:) 
            d(:,p)=d(:,p)-d2(:)-d4(:) 
             
        end if 
         
    end do 
     
    delta_t=cfl*vol/scale_a 
     
    return 
     
    end subroutine 
     
    subroutine r_smooth 
     
    use global 
    implicit none 
     
    integer::i 
     
    do i=1,ncells 
         
        r(:,i)=(d(:,i)-q(:,i))/vol(i) 
         
    end do 
     
    return 
     
    end subroutine 
     
    subroutine results 
     
    use global 
    implicit none 
     
    integer::i,j,k 
    real*8::u(nnodes),v(nnodes),rou(nnodes),p(nnodes),cp(nnodes) 
    real*8::mach(nnodes),nvol(nnodes),uv,ck,pk,fx,fy,cl,cd 
     
    u=0. 
    v=0. 
    rou=0. 
    p=0. 
    mach=0. 
    nvol=0. 
    fx=0. 
    fy=0. 
     
    do i=1,ncells 
         
        do j=1,nedges 
             
            if(iedge(3,j)==i .or. iedge(4,j)==i) then 
                 
                nvol(iedge(1,j))=nvol(iedge(1,j))+vol(i) 
                nvol(iedge(2,j))=nvol(iedge(2,j))+vol(i) 
                 
            end if 
             
        end do 
         
        do j=1,4 
             
            k=icell(j,i) 
            u(k)=u(k)+vol(i)*w(2,i)/w(1,i) 
            v(k)=v(k)+vol(i)*w(3,i)/w(1,i) 
            rou(k)=rou(k)+w(1,i)*vol(i) 
            pk=(gama-1.)*(w(4,i)-0.5*(w(2,i)**2+w(3,i)**2)/w(1,i)) 
            ck=sqrt(gama*pk/w(1,i)) 
            p(k)=p(k)+vol(i)*pk 
            mach(k)=mach(k)+vol(i)*sqrt((w(2,i)/w(1,i))**2+(w(3,i)/w(1,i))**2)/ck 
             
        end do 
         
    end do 
    nvol=nvol/2. 
    u=u/nvol 
    v=v/nvol 
    rou=rou/nvol 
    p=p/nvol 
    mach=mach/nvol 
    uv=0.5*gama*pinf*machinf**2 
    cp=(p-pinf)/uv 
     
    do i=1,nedges 
         
        if(iedge(4,i)==-1) then 
             
            fy=fy-0.5*(p(iedge(1,i))+p(iedge(2,i)))*dxy(1,i) 
            fx=fx+0.5*(p(iedge(1,i))+p(iedge(2,i)))*dxy(2,i) 
             
        end if 
         
    end do 
     
    cl=(fx*sin(angle)+fy*cos(angle))/uv 
    cd=(fx*cos(angle)-fy*sin(angle))/uv 
     
    open(30,file='results.dat') 
    write(30,*) 'TITLE="NACA0012 CFD RESULTS"' 
    write(30,*) 'VARIABLES="X","Y","U","V","<greek>r</greek>","P","MACH","C<sub>P</sub>"' 
    write(30,*) 'ZONE N=',nnodes,',E=',ncells,', F=FEPOINT, ET=quadrilateral' 
     
    do i=1,nnodes 
         
        write(30,*) xy(1,i),xy(2,i),u(i),v(i),rou(i),p(i),mach(i),cp(i) 
         
    end do 
     
    do i=1,ncells 
         
        write(30,*) icell(1,i),icell(2,i),icell(3,i),icell(4,i) 
         
    end do 
     
    close(30) 
     
    open(40,file='cl-cd.txt') 
    write(40,*) 'CL=',cl 
    write(40,*) 'CD=',cd 
    close(40) 
     
    open(50,file='cp.dat') 
     
    do i=1,nwing 
         
        write(50,*) xy(1,wing(1,i)),cp(wing(1,i)),cp(wing(2,i)) 
         
    end do 
     
    close(50) 
     
    return 
     
    end subroutine