如何采用过控制体积法来确定熔体前沿位置和形状的？

hyacinth · 发表于 2005-7-23 14:45

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

所指要的界面是指熔体与熔体、熔体与气体之间的界面，这些界面随着流体的流动而移动。我了解到有采用VOF(流体体积法）的，也有采用CV（控制体积法），CV是在VOF上发展起来的。不知是否有人做过这方面的编程？ 参考程序 8.1. 1D FORTRAN code. A sample code for solving 1D hyperbolic system of conservation laws is provided below. The users are asked to provide their own PDE solvers, such as ENO schemes, central schemes and BGK methods, to form a subroutine PDES . Inthe following code, GAMA is the ratio of specied heat capacities of the uid, GAM1 = GAMA-1 , and IADAP is a switch parameter. If IADAP < 0, then PDEs are solved on uniform mesh, otherwise on adaptive mesh. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C 1D FORTRAN Code of Adaptive Algorithm C Main Program CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PARAMETER(NX=3000) COMMON/COM1/UO(3,NX),UN(3,NX) & ,XO(NX),XN(NX) COMMON/COM2/NPX,NXI COMMON/EULER/GAMA,GAM1 DATA NPX,NXI/103,103/ DATA XLEFT,XRIGHT,XPO,GAMA & ,GAM1/0.,1.,.5,1.4,.4/ TSTOP=0.2 C C Initial partition and initial data C X: x-coordinate, old, N: new DX=(XRIGHT-XLEFT)/(NPX-3) DO J=1,NPX+4 XO(J)=XLEFT+(J-3)*DX XN(J)=XO(J) ENDDO CALL INITIAL_DATA WRITE(*,*)' IADAP<0 or >0' READ(*,*) IADAP C C Move mesh and solve PDEs IT=0 TIME=0.0 100 IF(TIME.GT.TSTOP) GOTO 200 CALL AMMM1D(IADAP) CALL PDES(DT) IT=IT+1 TIME=TIME+DT WRITE(*,*) IT,TIME GOTO 100 C Output the final results 200 CALL OUTPUT END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Mesh-moving & solution-updating CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE AMMM1D(IADAP) PARAMETER(NX=3000) IMPLICIT REAL*8 (A-H,O-Z) COMMON/COM1/UO(3,NX),UN(3,NX) & ,XO(NX),XN(NX) COMMON/COM2/NPX,NXI COMMON/EULER/GAMA,GAM1 DIMENSION WO(NNX), WN(NNX) &,DUX(3,NX),DUU(NX),DUS(NX) DATA ALPHA,BETA/100.,100./ FMON(X,DX,ALPHA,BETA)=SQRT(1. & +ALPHA*X*X+BETA*DX*DX) IF(IADAP.LT.0) RETURN DO J=1,NXI+2 XN(J)=XO(J) DO L=1,3 UN(L,J)=UO(L,J) ENDDO ENDDO IT=0 900 IT=IT+1 C C DEFINE MONITOR FUNCTION SUU=0.0 SUS=0.0 DO J=2,NXI+1 RHOP1=UO(1,J+1) UXP1 =UO(2,J+1)/RHOP1 PREP1=GAM1*(UO(3,J+1)-.5 & *RHOP1*UXP1*UXP1) RHOM1=UO(1,J-1) UXM1 =UO(2,J-1)/RHOM1 PREM1=GAM1*(UO(3,J-1)-0.5 & *RHOM1*UXM1*UXM1) DUU(J)=(UXP1-UXM1)/2.0 DUS(J)=((PREP1/RHOP1**GAMA) & -(PREM1/RHOM1**GAMA))/2.0 SUU=MAX(SUU,DUU(J)) SUS=MAX(SUS,DUS(J)) ENDDO DO J=2,NXI+1 UU=DUU(J)/(SUU+1.0E-10) DU=DUS(J)/(SUS+1.0E-10) WO(J)=FMON(UU,DU,ALPHA,BETA) ENDDO DO 100 IJK=1,4 DO J=2,NXI+1 WN(J)=WO(J) IF(J.GT.2.and.J.LT.NXI+1) THEN WN(J)=(WO(J+1)+2.*WO(J) & +WO(J-1))/4. ENDIF ENDDO DO J=3,NXI WO(J)=WN(J) ENDDO 100 CONTINUE C C update mesh DO J=3,NXI XN(J)=(WO(J)*XO(J+1)+WO(J-1) & *XN(J-1))/(WO(J)+WO(J-1)) ENDDO C C compute the iteration error SUM=0.0 DO J=1,NXI SUM=SUM+ABS(XN(J)-XO(J)) ENDDO SUM=SUM/NXI C C update solution using upwinding flux DO L=1,3 DO J=3,NXI DUL=UO(L,J)-UO(L,J-1) DUR=UO(L,J+1)-UO(L,J) DUX(L,J)=(SIGN(1.,DUL)+SIGN(1., & DUR))*ABS(DUL*DUR)/(ABS(DUL) & +ABS(DUR)+1.E-10) ENDDO DUX(L,2)=0.0 DUX(L,NXI+1)=0.0 ENDDO C C XIX: area of control volume C old; N: new DO J=NIA,NIB XIXO=XO(J+1)-XO(J) XIXN=XN(J+1)-XN(J) SP1=(XO(J+1)-XN(J+1))/XIXN SP2=(XO(J)-XN(J))/XIXN DO L=1,3 UJP1=(UO(L,J+1)-0.5*DUX(L,J+1)) UJP0=UO(L,J )-0.5*DUX(L,J) UJM1=UO(L,J-1)+0.5*DUX(L,J-1) UJM0=UO(L,J )+0.5*DUX(L,J) UN(L,J)=UO(L,J)*XIXO/XIXN-.5*( & SP1*(UJP1+UJM0)-SP2*(UJP0+UJM1)) & +.5*(ABS(SP1)*(UJP1-UJM0) & -ABS(SP2)*(UJP0-UJM1)) ENDDO ENDDO IF(SUM.GE.1.0E-6) THEN DO J=3,NXI DO L=1,3 UO(L,J)=UN(L,J) ENDDO XO(J)=XN(J) ENDDO GOTO 900 ENDIF RETURN END 8.2. 2D code. A sample code for solving 2D Riemann problem is provided below. Again the users are asked to provide their own PDE solver to form a PDE solution subroutine PDES . CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C 2D FORTRAN Code of Adaptive Algorithm C Main Program CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PARAMETER(NX=100,NY=100) COMMON/C1/UO(4,NX,NY),UN(4,NX,NY) COMMON/C2/XO(NX,NY),YO(NX,NY) COMMON/C3/NPX,NPY,NXI,NET COMMON/C4/GAMA,PI COMMON/C5/RHOA,UXA,UYA,PREA & ,RHOB,UXB,UYB,PREB DATA NPX,NPY,NXI,NET/53,53,53,53/ DATA XLEFT,YLEFT,XRIGHT,YRIGHT & ,GAMA/0.,0.,1.,1.,1.4/ PI=ASIN(1.0)*2.0 TSTOP=0.2 C C Initial partition and initial data DX=(XRIGHT-XLEFT)/(NPX-3) DY=(YRIGHT-YLEFT)/(NPY-3) DO J=1,NPX+4 DO K=1,NPY+4 XO(J,K)=XLEFT+(J-3)*DX YO(J,K)=YLEFT+(K-3)*DY XN(J,K)=XO(J,K) YN(J,K)=YO(J,K) ENDDO ENDDO CALL INITIAL_DATA WRITE(*,*) ' IADAP<0 or > 0' READ(*,*) IADAP C C Move mesh and solve PDEs IT=0 TIME=0.0 100 IF(TIME.GT.TSTOP) GOTO 200 CALL AMMM1D(IADAP) CALL PDES(DT) IT=IT+1 TIME=TIME+DT WRITE(*,*) IT,TIME GOTO 100 C C Output the final results 200 CALL OUTPUT END SUBROUTINE AMMM2D(IADAP) PARAMETER(NX=100,NY=100) COMMON/C1/UO(4,NX,NY),UN(4,NX,NY) COMMON/C2/XO(NX,NY),YO(NX,NY) COMMON/C3/NPX,NPY,NXI,NET COMMON/C4/GAMA,PI COMMON/C5/RHOA,UXA,UYA,PREA & ,RHOB,UXB,UYB,PREB DIMENSION WO(NX,NY),WN(NX,NY), & FLUX(4,NX,NY),FLUY(4,NX,NY) DIMENSION XN(NX,NY),YN(NX,NY), & DUX(4,NNX) DATA ALPHA,BETA,EP/0.,500.,1.E-10/ FMON(x,DX,DY,ALPHA,BETA)=SQRT(1.0 & +ALPHA*X*X+BETA*(DX*DX+DY*DY) ) IF(IADAP.LT.0) RETURN NIA=4 NIB=NXI-1 NTA=4 NTB=NET-1 NXA=3 NXB=NPX-1 NYA=3 NYB=NPY-1 DO 100 J=1,NXI+2 DO 100 K=1,NET+2 XN(J,K)=XO(J,K) YN(J,K)=YO(J,K) 100 CONTINUE DO 110 J=NXA-2,NXB+2 DO 110 K=NYA-2,NYB+2 DO 110 L=1,4 UN(L,J,K)=UO(L,J,K) 110 CONTINUE IT=0 120 IT=IT+1 C C give monitor function DO 130 J=NXA-1,NXB+1 DO 130 K=NYA-1,NYB+1 UU=UO(1,J,K) DUXI=.5*(UO(1,J+1,K)-UO(1,J-1,K)) DUET=.5*(UO(1,J,K+1)-UO(1,J,K-1)) WO(J,K)=FMON(UU,DUXI,DUET & ,ALPHA,BETA) 130 CONTINUE DO 160 IJK=1,4 DO 140 J=NXA-1,NXB+1 DO 140 K=NYA-1,NYB+1 WN(J,K)=WO(J,K) IF(J.GE.NIA.and.J.LE.NIB) THEN IF(K.GE.NTA.and.K.LE.NTB) THEN WN(J,K)=WO(J,K)/4.+(WO(J+1,K) & +WO(J-1,K)+WO(J,K+1)+WO(J,K-1))/8. & +(WO(J+1,K+1)+WO(J+1,K-1) & +WO(J-1,K+1)+WO(J-1,K-1))/16. ENDIF ENDIF 140 CONTINUE DO 150 J=NXA-1,NXB+1 DO 150 K=NYA-1,NYB+1 WO(J,K)=WN(J,K) 150 CONTINUE 160 CONTINUE DO 170 J=NIA,NIB DO 170 K=NTA,NTB CXP=0.5*(WO(J,K)+WO(J,K-1)) CXM=0.5*(WO(J-1,K)+WO(J-1,K-1)) CYP=0.5*(WO(J,K)+WO(J-1,K)) CYM=0.5*(WO(J,K-1)+WO(J-1,K-1)) CXY=(CXP+CXM+CYP+CYM) XN(J,K)=( CXP*XO(J+1,K)+CXM*XN(J-1,K) & +CYP*XO(J,K+1)+CYM*XN(J,K-1) )/CXY YN(J,K)=( CXP*YO(J+1,K)+CXM*YN(J-1,K) & +CYP*YO(J,K+1)+CYM*YN(J,K-1) )/CXY 170 CONTINUE C C boundary grid redistribution DO K=NTA,NTB J=NIA-1 XN(J,K)=XO(J,K) YN(J,K)=YO(J,K)+(YN(J+1,K)-YO(J+1,K)) J=NIB+1 XN(J,K)=XO(J,K) YN(J,K)=YO(J,K)+YN(J-1,K)-YO(J-1,K) ENDDO DO J=NIA,NIB K=NTA-1 YN(J,K)=YO(J,K) XN(J,K)=XO(J,K)+XN(J,K+1)-XO(J,K+1) K=NTB+1 YN(J,K)=YO(J,K) XN(J,K)=XO(J,K)+XN(J,K-1)-XO(J,K-1) ENDDO C SUM=0.0 DO 180 J=NIA,NIB DO 180 K=NTA,NTB SUM=SUM+ABS(XO(J,K)-XN(J,K)) & +ABS(YO(J,K)-YN(J,K)) 180 CONTINUE SUM=SUM/(NIB-NIA)/(NTB-NTA) C C x-direction flux DO 190 K=NYA,NYB DO L=1,4 DO J=NXA,NXB DUL=(UO(L,J,K)-UO(L,J-1,K)) & /(XO(J,K)-XO(J-1,K)) DUR=(UO(L,J+1,K)-UO(L,J,K)) & /(XO(J+1,K)-XO(J,K)) DUX(L,J)=(SIGN(1.,DUL) & +SIGN(1.,DUR))*ABS(DUL & *DUR)/(ABS(DUL)+ABS(DUR)+EP) ENDDO DUX(L,NXA-1)=0.0 DUX(L,NXB+1)=0.0 ENDDO DO J=NXA-1,NXB DTX=XO(j+1,k+1)-XO(j+1,k ) DTY=YO(j+1,k+1)-YO(j+1,k ) ALENG=SQRT(DTX**2+DTY**2) SZTAK=DTY/ALENG CZTAK=DTX/ALENG SPX=.5*(XO(J+1,K)+XO(J+1,K+1) & -XN(J+1,K)-XN(J+1,K+1)) SPY=.5*(YO(J+1,K)+YO(J+1,K+1) & -YN(J+1,K)-YN(J+1,K+1)) SPN=SZTAK*SPX-CZTAK*SPY DO L=1,4 UUL=UO(L,J,K)+.5*(XO(J+1,K) & -XO(J,K))*DUX(L,J) UUR=UO(L,J+1,K)-.5*(XO(J+1,K) & -XO(J,K))*DUX(L,J+1) FLUX(L,J,K)=.5*(SPN*(UUR+UUL) & -ABS(SPN)*(UUR-UUL))*ALENG ENDDO ENDDO 190 CONTINUE C C y-direction flux DO 200 J=NXA,NXB DO L=1,4 DO K=NYA,NYB DUL=(UO(L,J,K)-UO(L,J,K-1)) & /(YO(J,K)-YO(J,K-1)) DUR=(UO(L,J,K+1)-UO(L,J,K)) & /(YO(J,K+1)-YO(J,K)) DUX(L,K)=(SIGN(1.,DUL)+ & SIGN(1.,DUR))*ABS(DUL*DUR) & /(ABS(DUL)+ABS(DUR)+EP) ENDDO DUX(L,NYA-1)=0.0 DUX(L,NYB+1)=0.0 ENDDO DO K=NYA-1,NYB DTX=XO(j,k+1)-XO(j+1,k+1) DTY=YO(j,k+1)-YO(j+1,k+1) ALENG=SQRT(DTX**2+DTY**2) SZTAK=DTY/ALENG CZTAK=DTX/ALENG SPX=.5*(XO(J,K+1)+XO(J+1,K+1) & -XN(J,K+1)-XN(J+1,K+1)) SPY=.5*(YO(J,K+1)+YO(J+1,K+1) & -YN(J,K+1)-YN(J+1,K+1)) SPN=SZTAK*SPX-CZTAK*SPY DO L=1,4 UUL=UO(L,J,K)+.5*(YO(J,K+1) & -YO(J,K))*DUX(L,K) UUR=UO(L,J,K+1)-.5*(YO(J,K+1) & -YO(J,K))*DUX(L,K+1) FLUY(L,J,K)=.5*(SPN*(UUR+UUL) & -ABS(SPN)*(UUR-UUL))*ALENG ENDDO ENDDO 200 CONTINUE C C areas of control volumes C and update the solution values DO 210 K=NYA,NYB DO 210 J=NXA,NXB AREAO=.5*((XO(j,k)-XO(j+1,k+1)) & *(YO(j+1,k)-YO(j,k+1)) & -(YO(j,k)-YO(j+1,k+1)) & *(XO(j+1,k)-XO(j,k+1))) AREAN=0.5*((XN(j,k)-XN(j+1,k+1)) & *(YN(j+1,k)-YN(j,k+1)) & - (YN(j,k)-YN(j+1,k+1)) & *(XN(j+1,k)-XN(j,k+1))) DO L=1,4 UNN=UO(L,J,K)*AREAO-(FLUX(L,J,K) & -FLUX(L,J-1,K))-(FLUY(L,J,K) & -FLUY(L,J,K-1)) UN(L,J,K)=UNN/AREAN ENDDO 210 CONTINUE IF (SUM.GE.1.E-5) THEN DO 220 J=NXA,NXB DO 220 K=NYA,NYB DO 220 L=1,4 UO(L,J,K)=UN(L,J,K) 220 CONTINUE DO 230 J=NXA-1,NXB+1 DO 230 K=NYA-1,NYB+1 XO(J,K)=XN(J,K) YO(J,K)=YN(J,K) 230 CONTINUE GOTO 120 ENDIF RETURN END

FSI · 发表于 2005-7-24 13:31

这属于界面追踪方法吧，有用 LevelSet 作的。
你贴的代码是不是根据计算结果来移动网格点？

[此贴子已经被作者于2005-7-24 13:33:26编辑过]

ik760218 · 发表于 2005-7-26 00:26

曾经看过FE/CV法的基本思路，它是使用有限元求解压力场，而控制体积法用于确定熔体前沿位置和形状。具体步骤是首先将模腔充填区域用有限元技术划分为有限单元，而在有限单元的各节点上再构造附加的控制体积元，然后将整个充模瞬态过程离散为一系列亚稳态过程。求解压力场，再由压力场计算速度场，从速度场计算流量，最后根据Courant条件估算最小时间步长，并计算用于追踪熔体前沿的填充因子。整个计算过程不断重复直到填满模腔

账号		自动登录	找回密码
密码			我要加入

[基础理论] 如何采用过控制体积法来确定熔体前沿位置和形状的？

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。