I played a bit around with the laws of physics, and when setting up a system with small light and big heavy balls/particles, I managed to get something that looks like a liquid and a gas phase. It can be adjusted in quite a lot of ways. You are welcome to make a use of it if you can find one, play around with it a bit if you like.
SeedRnd(1) Graphics 1024,768 SetBuffer BackBuffer() Const ballno=199 ;Const gravg#=0.0001 Dim x#(ballno) Dim y#(ballno) Dim dx#(ballno) Dim dy#(ballno) Dim mass#(ballno) Dim radius#(ballno) For a=0 To ballno x(a)=Rand(512-50-50+25,512+50+50-25) y(a)=Rand(384-50-50+25,384+50+50-25)+250 ; dx(a)=Rnd(-1,1) ; dy(a)=Rnd(-1,1) mass(a)=(a/40)^3+1 radius(a)=mass(a)^(1/3.0)*1 Next ;mass(0)=1 ;x(0)=200 ;dx(0)=.5 ;y(0)=300 ;x(1)=600 ;dx(1)=-.5 ;y(1)=300 ;mass(1)=1 While Not KeyHit(1) ; For c=0 To 3 For a=0 To ballno x(a)=x(a)+dx(a) y(a)=y(a)+dy(a)+.0005 dy(a)=dy(a)+.001 If x(a)<10+400 And dx(a)<0 Then dx(a)=-dx(a) ; If y(a)<10 And dy(a)<0 Then dy(a)=-dy(a) If x(a)>1013-400 And dx(a)>0 Then dx(a)=-dx(a) If y(a)>757-50 And dy(a)>0 Then dy(a)=-dy(a) Next For a=0 To ballno-1 For b=a+1 To ballno dist#=(x(a)-x(b))^2+(y(a)-y(b))^2 sqrdist#=Sqr(dist)*1000 dx(a)=dx(a)+(x(b)-x(a))*mass(b)/dist/sqrdist dy(a)=dy(a)+(y(b)-y(a))*mass(b)/dist/sqrdist dx(b)=dx(b)+(x(a)-x(b))*mass(a)/dist/sqrdist dy(b)=dy(b)+(y(a)-y(b))*mass(a)/dist/sqrdist distmin=(radius(a)+radius(b))^2 If dist<distmin Then If (x(a)+dx(a)-x(b)-dx(b))^2+(y(a)+dy(a)-y(b)-dy(b))^2<dist Then atobx#=x(b)-x(a) atoby#=y(b)-y(a) lenatob2#=atobx^2+atoby^2 damdbx#=dx(a)-dx(b) damdby#=dy(a)-dy(b) factora#=-(damdbx*atobx+damdby*atoby)/lenatob2 factorb#=(damdbx*atobx+damdby*atoby)/lenatob2 aonatobx#=atobx*factora aonatoby#=atoby*factora bonatobx#=atobx*factorb bonatoby#=atoby*factorb dx(a)=dx(a)+aonatobx*mass(b)/(mass(a)+mass(b))*2 dy(a)=dy(a)+aonatoby*mass(b)/(mass(a)+mass(b))*2 dx(b)=dx(b)+bonatobx*mass(a)/(mass(a)+mass(b))*2 dy(b)=dy(b)+bonatoby*mass(a)/(mass(a)+mass(b))*2 ; While (x(a)+dx(a)-x(b)-dx(b))^2+(y(a)+dy(a)-y(b)-dy(b))^2<dist ; x(a)=x(a)+dx(a)/2 ; y(a)=y(a)+dy(a)/2+.0025 ; dy(a)=dy(a)+.005 ; x(b)=x(b)+dx(b)/2 ; y(b)=y(b)+dy(b)/2+.0025 ; dy(b)=dy(b)+.005 ; Wend End If End If dx(a)=dx(a)+(x(b)-x(a))*mass(b)/dist/sqrdist dy(a)=dy(a)+(y(b)-y(a))*mass(b)/dist/sqrdist dx(b)=dx(b)+(x(a)-x(b))*mass(a)/dist/sqrdist dy(b)=dy(b)+(y(a)-y(b))*mass(a)/dist/sqrdist Next Next ; Next Cls For a=0 To ballno Oval x(a)-radius(a),y(a)-radius(a),radius(a)*2,radius(a)*2 Next Flip Wend