求救壁面函数！哪位好心老师能帮帮我

hhf-dxj

我是长沙理工大学的一位研究生，请各位老师帮忙，请问各位老师有没有壁面函数的详细的推导过程和其源程序，我在做毕业论文，真的很需要用，但是我不会，麻烦各位了，如果哪位好心的老师有麻烦将其传上来，晚辈在此表示深深的谢意了！

harry

      SUBROUTINE BSWALL(IFLAG)                                                  
C     SOLID WALL BOUNDARIES:                                                   
C     IFLAG = 1(NORTH), 2(SOUTH), 3(EAST), 4(WEST)           
C***********************************************************************        
      INCLUDE ';COM3D1.FI';                                                         
C                                                                              
        MFLAG1=0                                                               
        MFLAG2=0                                                               
        MFLAG3=0                                                               
        MFLAG4=0                                                               
      IF(     IFLAG.EQ.1) THEN                                                  
        MFLAG1=1                                                               
      ELSE IF(IFLAG.EQ.2) THEN                                                  
        MFLAG2=1                                                               
      ELSE IF(IFLAG.EQ.3) THEN                                                  
        MFLAG3=1                                                               
      ELSE IF(IFLAG.EQ.4) THEN                                                  
        MFLAG4=1                                                               
      ENDIF                                                                     
C                                                                              
        MWAL = MWALN*MFLAG1 + MWALS*MFLAG2 + MWALE*MFLAG3                       
     >       + MWALW*MFLAG4                     
      ISANOW =  ISAN*MFLAG1 +  ISAS*MFLAG2 +  ISAE*MFLAG3                       
     >       +  ISAW*MFLAG4                       
      DO 1110 IM=1,MWAL                                                         
          IIP= INDWN(IM)*MFLAG1 + INDWS(IM)*MFLAG2 + INDWE(IM)*MFLAG3           
     >       + INDWW(IM)*MFLAG4            
      IIB=IIP+        NI*MFLAG1 -        NI*MFLAG2 +           MFLAG3           
     >       -           MFLAG4           
      ID1=IIP+    (1+NI)*MFLAG1                    +           MFLAG3           
     >       +        NI*MFLAG4            
      ID2=IIP+        NI*MFLAG1 +           MFLAG2 +    (NI+1)*MFLAG3           
C     NORMAL VECTOR, AREA AND NORMAL DISTANCE                                   
             DAX=F(ISCOX+ID1)-F(ISCOX+ID2)                                      
             DAY=F(ISCOY+ID1)-F(ISCOY+ID2)                                      
         DN1=-DAY                                                   
         DN2= DAX                                                   
        AREA=0.5*SQRT(DN1*DN1+DN2*DN2)                                 
          A1=0.5*DN1/(AREA+SMALL)                                               
          A2=0.5*DN2/(AREA+SMALL)                                               
       DELTA=ABS((F(ISX+IIP)-F(ISX+IIB))*A1+(F(ISY+IIP)-F(ISY+IIB))*A2      
     >          )+SMALL                              
C     SHEAR STRESS AT THE WALL                                                  
            A11=1.-A1*A1                                                        
            A22=1.-A2*A2                                                        
            A12=A1*A2                                                           
         VXCOMP=A11*F(ISU+IIP)-A12*F(ISV+IIP)                    
         VYCOMP=A22*F(ISV+IIP)-A12*F(ISU+IIP)                    
          VPARL=SQRT(VXCOMP*VXCOMP+VYCOMP*VYCOMP)                 
          SQRTK=SQRT(ABS(F(ISTE+IIP)))                                          
          YPLUS=F(ISDEN+IIP)*SQRTK*CDQR*DELTA/VISCOS                           
          TMULT=VISCOS/DELTA                                                   
         IF(YPLUS.GT.11.6)                                                      
     >    TMULT=F(ISDEN+IIP)*CDQR*SQRTK*CAPPA/ALOG(ELOG*YPLUS)                  
C     MODIFICATION OF SOURCE TERMS                                             
        FACTOR=TMULT*AREA                                                      
          UPUB=F(ISU+IIP)-F(ISU+IIB)                                            
          VPVB=F(ISV+IIP)-F(ISV+IIB)                                            
      IF(    NN.EQ.2) THEN                                                      
       F(ISSP+IIP)=F(ISSP+IIP)-FACTOR*A11                                       
       F(ISSU+IIP)=F(ISSU+IIP)+FACTOR*(F(ISU+IIB)*A11+VPVB*A12)        
      ELSEIF(NN.EQ.3) THEN                                                      
       F(ISSP+IIP)=F(ISSP+IIP)-FACTOR*A22                                       
       F(ISSU+IIP)=F(ISSU+IIP)+FACTOR*(F(ISV+IIB)*A22+UPUB*A12)        
      ELSEIF(NN.EQ.4) THEN                                                      
          IF(YPLUS.LE.11.6) THEN                                                
            TERM=F(ISDEN+IIP)*CDTQ*SQRTK*YPLUS/DELTA                           
           ELSE                                                                 
            TERM=F(ISDEN+IIP)*CDTQ*SQRTK*ALOG(ELOG*YPLUS)/(CAPPA*DELTA)         
          ENDIF                                                                 
          TAUW=TMULT*VPARL                                                      
          GENR=TAUW*TAUW/(CAPPA*VISCOS*YPLUS+SMALL)                             
       F(ISSP+IIP)=TESPC(IIP)-TERM*F(ISVOLP+IIP)                                
       F(ISSU+IIP)=TESUC(IIP)+GENR*F(ISVOLP+IIP)                                
      ELSEIF(NN.EQ.5) THEN                                                      
       F(ISSP+IIP)=-GREAT                                                      
       F(ISSU+IIP)= GREAT*ABS(F(ISTE+IIP))**1.5*CDTQ/(CAPPA*DELTA)              
CTANG.TANG.TANG.--------------------------------------------------
      ELSEIF(NN.EQ.6) THEN
       F(ISS+IIB)=F(ISS+IIP)
CTANG.TANG.TANG.--------------------------------------------------
      ENDIF                                                                     
1110 F(ISANOW+IIP)=0.0                                                         
      RETURN                                                                    

harry

你可以参考peric的程序，上边这个是实施坐标变换的2d程序（peric的程序没实施坐标变换，直接在物理网格上求解），peric的书上推得很详细