粒子法で、ざっぱ〜ん！

時間間隔がまずかったみたいで、ちゃんとざっぱ〜んとなりました。
最後爆発しちゃいますが。

あと、実際には粒子の速さによって計算する時間間隔を調整しないといけなかったり、粒子が近づき過ぎないようにする処理が必要になるのですが、まぁ適当に動いたのでとりあえずここまで。これらの処理を入れると爆発しなくなるんじゃないかと思います。
Flashはこっち。
http://rps13.180.jp/parc3.swf

粒子法でやってます。このPDFに解説が。
http://www.nagare.or.jp/nagare/21-3/21-3-t02.pdf
概念はこのプレゼンがわかりやすい
http://www.prometech.co.jp/uploads/img/MPS_koushukai_shibata.pdf
ここにムービーがあります。「ここでは2つほど紹介します。」のリンク
http://www.geolab.jp/Feel&Think/2005/FT_sato0505.html

参考にした本はこれです。

粒子法 (計算力学レクチャーシリーズ)

ほぼ、本の数式どおりですが、いくつか、付属CDのサンプルソースから処理を持ってきました。ここは本にも解説がみあたらず、何してるかわかりません。

この本には粒子法のMPEGムービーがついてます。

パソコンで見る流れの科学―数値流体力学入門 CD-ROM付 (ブルーバックス)

ソースはこんな感じで。ちょっと長いですが。

import java.awt.*;
import java.awt.image.BufferedImage;
import javax.swing.*;

public class Particle {
    enum Type{
        WATER(Color.BLUE), OUTWALL(Color.BLACK), INWALL(Color.GREEN);

        Color color;
        Type(Color c){
            color = c;
        }
    }
    
    static class Elem{
        double x;
        double y;
        double vx;
        double vy;
        double p;
        Type type;
        boolean surface = false;

        public Elem(double x, double y, Type type) {
            this.x = x;
            this.y = y;
            this.type = type;
        }
        Color getColor(){
            return surface ? Color.WHITE : type.color;
        }
    }
    
    static Elem[] elements = new Elem[1000];
    static int elementsCount;
    static JLabel label;
    static Image img = new BufferedImage(500, 400, BufferedImage.TYPE_INT_RGB);
    
    static final double parcSize = 1.5 / 100;//粒子サイズ
    static final double re = 2.1;//近傍粒子の範囲
    static final double relap = 4.0;//ラプラシアン用の近傍粒子の範囲
    static final double deltat = 0.0001;//時間の単位
    static final double grav = 9.8;//重力加速度(m/s^2)
    static final double dim = 2;//次元数
    static final double nyu = 1.0038E-6;//動粘性係数(水 1.0038E-6 m^2/s)
    static final double row = 1.0e3;

    static double n0;//基準粒子密度
    static double lambda;//ラプラシアン用の係数
    static final double beta = 0.95;//自由表面条件
    
    public static void main(String[] args){
        //左右の外壁
        double maxy = 0;
        for(double y = 0; y < 0.5; y += parcSize){
            elements[elementsCount++] = new Elem(0, y, Type.INWALL);
            elements[elementsCount++] = new Elem(-parcSize, y, Type.OUTWALL);
            elements[elementsCount++] = new Elem(-parcSize * 2, y, Type.OUTWALL);
            elements[elementsCount++] = new Elem(.8, y, Type.INWALL);
            elements[elementsCount++] = new Elem(.8+parcSize, y, Type.OUTWALL);
            elements[elementsCount++] = new Elem(.8+parcSize * 2, y, Type.OUTWALL);
            maxy = y;
        }
        //下の壁
        for(double x = -parcSize * 2; x < .8 + parcSize * 2; x += parcSize){
            elements[elementsCount++] = new Elem(x, maxy + parcSize, 
                    (x < 0 || x > .8) ? Type.OUTWALL : Type.INWALL);
            elements[elementsCount++] = new Elem(x, maxy + parcSize * 2, Type.OUTWALL);
            elements[elementsCount++] = new Elem(x, maxy + parcSize * 3, Type.OUTWALL);
        }
        
        //水柱
        int sample = 0;
        for(double x = parcSize; x < .30; x+= parcSize){
            for(double y = parcSize * 3; y < .50; y += parcSize){
                if(y > .20 && y < .25 && x > .10 && x < .15){
                    //内部になる点を一点記憶しておく
                    sample = elementsCount;
                }
                elements[elementsCount++] = new Elem(x, y, Type.WATER);
            }
        }

        //基本の粒子密度n0とラプラシアン用の係数rambda
        n0 = 0;
        double a = 0;
        double b = 0;
        for(int i = 0; i < elementsCount; ++i){
            if(i == sample) continue;

            double d = Math.sqrt(
                    (elements[sample].x - elements[i].x) * (elements[sample].x - elements[i].x) +
                    (elements[sample].y - elements[i].y) * (elements[sample].y - elements[i].y));
            double w = (d < re * parcSize) ? parcSize * re / d - 1 : 0;            
            n0 += w;
            a += (elements[sample].x - elements[i].x) * (elements[sample].x - elements[i].x) +
                    (elements[sample].y - elements[i].y) * (elements[sample].y - elements[i].y) * w;
        }
        lambda = a / n0;

        //ウィンドウの表示
        JFrame f = new JFrame("粒子法");
        f.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
        f.setSize(500, 450);
        label = new JLabel(new ImageIcon(img));
        draw();
        f.add(label);
        f.setVisible(true);

        //処理スレッド
        new Thread(){
            @Override
            public void run() {
                for(;;){
                    calc();
                    draw();
                }
            }
        }.start();
    }
    
    static void calc(){
        double[][] w = new double[elementsCount][elementsCount];
        //粒子密度を計算しておく
        for(int i = 0; i < elementsCount; ++i){
            for(int j = 0; j < i; ++j){
                double d = Math.sqrt(
                        (elements[i].x - elements[j].x) * (elements[i].x - elements[j].x) +
                        (elements[i].y - elements[j].y) * (elements[i].y - elements[j].y));
                double wd = (d < relap * parcSize) ? parcSize * relap / d - 1 : 0;
                w[i][j] = wd;
                w[j][i] = wd;
            }
        }
        
        //仮の速度・位置を計算
        for(int i = 0; i < elementsCount; ++i){
            if(elements[i].type != Type.WATER){
                continue;
            }
                
            //速度のラプラシアンを計算
            double laplasux = 0;
            double laplasuy = 0;
            for(int j = 0; j < elementsCount; ++j){
                if(j == i) continue;
                if(w[i][j] == 0) continue;
                laplasux += (elements[j].vx - elements[i].vx) * w[i][j];
                laplasuy += (elements[j].vy - elements[i].vy) * w[i][j];
            }
            laplasux *= 2 * dim / n0 / lambda;
            laplasuy *= 2 * dim / n0 / lambda;
            
            //仮の位置・速度を計算
            elements[i].vx += + deltat * (nyu * laplasux);
            elements[i].vy += deltat * (nyu * laplasuy + grav);
            elements[i].x += elements[i].vx * deltat;
            elements[i].y += elements[i].vy * deltat;
        }

        //仮の粒子密度を計算
        double[][] wasta = new double[elementsCount][elementsCount];
        double[][] rsqu = new double[elementsCount][elementsCount];
        int[] neighborcount = new int[elementsCount];
        int[][] neighbor = new int[elementsCount][elementsCount];
        for(int i = 0; i < elementsCount; ++i){
            for(int j = 0; j < i; ++j){
                rsqu[i][j] = rsqu[j][i] =
                        (elements[i].x - elements[j].x) * (elements[i].x - elements[j].x) +
                        (elements[i].y - elements[j].y) * (elements[i].y - elements[j].y);
                double d = Math.sqrt(rsqu[i][j]);
                double rp = re * parcSize;
                if(d < rp){
                    double wd = rp / d - 1;
                    wasta[i][j] = wd;
                    wasta[j][i] = wd;
                    neighbor[i][neighborcount[i]++] = j;
                    neighbor[j][neighborcount[j]++] = i;
                }
            }
        }
        
        //粒子密度から自由表面になるところを検出
        double[] nasta = new double[elementsCount];
        for(int i = 0; i < elementsCount; ++i){
            for(int j = 0; j < neighborcount[i]; ++j){
                nasta[i] += wasta[i][neighbor[i][j]];
            }
            if(elements[i].type == Type.WATER){
                elements[i].surface = nasta[i] < n0 * beta;
            }
        }
        
        //よくわからんけどサンプルにあった
        int[] check = new int[elementsCount];
        for(int i = 0; i < elementsCount; ++i){
            if(elements[i].type != Type.WATER){
                check[i] = 2;
            }else if(elements[i].surface){
                check[i] = 1;
            }else{
                check[i] = 0;
            }
        }
        for(boolean cont = true; !cont;){
            cont = false;
            for(int i = 0; i < elementsCount; ++i){
                if(check[i] == 1){
                    for(int j = 0; j < neighborcount[i]; ++j){
                        if(check[neighbor[i][j]] == 0) check[neighbor[i][j]] = 1;
                    }
                    check[i] = 2;
                    cont = true;
                }
            }
        }
        
        //連立方程式の係数行列を作成
        double[] rk = new double[elementsCount];
        double[][] coef = new double[elementsCount][elementsCount];
        double[] val = new double[elementsCount];
        for(int i = 0; i < elementsCount; ++i){
            if(elements[i].surface){
                val[i] = 0;
            }else{
                val[i] = -row / deltat / deltat * (nasta[i] - n0) * lambda / 2 / dim;
            }
            for(int j = 0; j < neighborcount[i]; ++j){
                coef[i][j + 1] = wasta[i][neighbor[i][j]];
                coef[i][0] -= wasta[i][neighbor[i][j]];
            }
            coef[i][0] -= 1.0e-7 / deltat / deltat * row * n0;//よくわからないけどソースから
            if(check[i] == 0) coef[i][0] *= 2;
        }
        
        //連立方程式を解く(とりあえずガウスザイデル法。CG法がいいらしい)
        for(int n = 0; n < 1000; ++n){
            double div = 0;
            for(int i = 0; i < elementsCount; ++i){
                double sum = 0;
                for(int j = 0; j < neighborcount[i]; ++j){
                    sum += coef[i][j + 1] * rk[neighbor[i][j]];
                }
                
                double old = rk[i];
                rk[i] = (val[i] - sum) / coef[i][0];
                
                div += Math.abs(rk[i] - old);
            }
            //収束判定
            if(div < 0.0001) break;
        }
        
        
        for(int i = 0; i < elementsCount; ++i){
            //圧力が負の場合は0に
            if(rk[i] < 0) rk[i] = 0;
        }
        
        //速度・場所を再計算
        for(int i = 0; i < elementsCount; ++i){
            if(elements[i].type != Type.WATER) continue;//場所がかわるのは水だけ
            
            double udashx = 0;
            double udashy = 0;
            //圧力の最小値を得る
            double hatp = rk[i];
            for(int j = 0; j < neighborcount[i]; ++j){
                if(hatp > rk[neighbor[i][j]]) hatp = rk[neighbor[i][j]];
            }
            //速度の修正量を計算
            for(int j = 0; j < neighborcount[i]; ++j){
                int n = neighbor[i][j];
                udashx += (rk[n] - hatp) / rsqu[i][n] * (elements[n].x - elements[i].x) * wasta[i][n];
                udashy += (rk[n] - hatp) / rsqu[i][n] * (elements[n].y - elements[i].y) * wasta[i][n];
            }
            udashx = -udashx * deltat / row * dim / n0;
            udashy = -udashy * deltat / row * dim / n0;
            //速度の再計算
            elements[i].vx += udashx;
            elements[i].vy += udashy;
            
            //場所を再計算
            elements[i].x += deltat * udashx;
            elements[i].y += deltat * udashy;
        }
    }
    
    static void draw(){
        Graphics2D g = (Graphics2D) img.getGraphics();
        g.setColor(new Color(128, 212, 255));
        g.fillRect(0, 0, 500, 400);
        int size = 8;
        double s = size / parcSize;
        for(int i = 0; i < elementsCount; ++i){
            Elem el = elements[i];
            g.setColor(el.getColor());
            g.fillOval((int)(el.x * s) + 30, (int)(el.y * s) + 20, size, size);
        }
        g.dispose();
        label.repaint();
    }
}