声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 2752|回复: 2

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

[复制链接]
发表于 2005-7-23 14:45 | 显示全部楼层 |阅读模式

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

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

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

使用道具 举报

发表于 2005-7-24 13:31 | 显示全部楼层
<P>这属于界面追踪方法吧,有用 LevelSet 作的。</P>
<P>你贴的代码是不是根据计算结果来移动网格点?</P>
[此贴子已经被作者于2005-7-24 13:33:26编辑过]

发表于 2005-7-26 00:26 | 显示全部楼层
曾经看过FE/CV法的基本思路,它是使用有限元求解压力场,而控制体积法用于确定熔体前沿位置和形状。具体步骤是首先将模腔充填区域用有限元技术划分为有限单元,而在有限单元的各节点上再构造附加的控制体积元,然后将整个充模瞬态过程离散为一系列亚稳态过程。求解压力场,再由压力场计算速度场,从速度场计算流量,最后根据Courant条件估算最小时间步长,并计算用于追踪熔体前沿的填充因子。整个计算过程不断重复直到填满模腔
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-3-28 18:18 , Processed in 0.183166 second(s), 17 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表