简介

其实流体力学就是专有名词多, 比如流固耦合. 其实就是流体对固体的一个作用. 固体对流体的一个作用而已. 搞得神神秘秘的服了.

使用诺依曼边界条件的建议水滴代码

using UnityEngine;
using System.Collections;

public class wave_motion : MonoBehaviour 
{
	int size 		= 100;
	float rate 		= 0.005f; // like alpha ?
	float gamma		= 0.004f;
	float damping 	= 0.996f;
	float[,] 	old_h;
	float[,]	low_h; // ei
	float[,]	vh;
	float[,]	b;

	bool [,]	cg_mask;
	float[,]	cg_p;
	float[,]	cg_r;
	float[,]	cg_Ap;
	bool 	tag=true;

	Vector3 	cube_v = Vector3.zero;
	Vector3 	cube_w = Vector3.zero;

	int[,] dirs;
	// Use this for initialization
	void Start () 
	{
		Mesh mesh = GetComponent ().mesh;
		mesh.Clear ();

		Vector3[] X=new Vector3[size*size];

		for (int i=0; i=0 && j>=0 && i=0 && j>=0 && i=0 && j>=0 && i=0 && j>=0 && i=0 && j>=0 && i().mesh;
		Vector3[] X = mesh.vertices;

		int li = (int)Mathf.Ceil((p.x - 0.5f + 5f) / 0.1f);
		int lj = (int)Mathf.Ceil((p.z - 0.5f + 5f) / 0.1f);
		int ui = (int)Mathf.Floor((p.x + 0.5f + 5f) / 0.1f);
		int uj = (int)Mathf.Floor((p.z + 0.5f + 5f) / 0.1f);

		// for block 1, calculate low_h;
		for (int i = li; i <= ui; i++)
		for (int j = lj; j <= uj; j++) 
		if(i >= 0 && j >= 0 && i < size && j < size)
        {
			cg_mask[i, j] = true;
			low_h[i, j] = 0f;
        }
        else
        {
			cg_mask[i, j] = false;
			low_h[i, j] = new_h[i, j];
        }

		// then set up b and cg_mask for conjugate gradient.
		for(int i=0; i 0 ? h[i - 1, j] : h[i, j];
				float hi1j = i + 1 < size ? h[i + 1, j] : h[i, j];
				float hij_1 = j > 0 ? h[i, j - 1] : h[i, j];
				float hij1 = j + 1 < size ? h[i, j + 1] : h[i, j];
				new_h[i, j] = h[i, j] + (h[i, j] - old_h[i, j]) * damping + 
					rate * (hi1j + hij1 + hi_1j + hij_1 - 4 * h[i, j]);
			}
		}
		//Step 2: Block->Water coupling
		//for (int i = 0; i < size; i++)
		//{
		//	for (int j = 0; j < size; j++)
		//	{
		//		vh[i, j] = 0;
		//	}
		//}
		//TODO: for block 1, calculate low_h.
		//TODO: then set up b and cg_mask for conjugate gradient.
		//TODO: Solve the Poisson equation to obtain vh (virtual height).
		//get_vh(new_h, 1);
		//TODO: for block 2, calculate low_h.
		//TODO: then set up b and cg_mask for conjugate gradient.
		//TODO: Solve the Poisson equation to obtain vh (virtual height).
		//get_vh(new_h, 2);
		//TODO: Diminish vh.

		//TODO: Update new_h by vh.
		//for (int i = 0; i < size; i++)
		//	for (int j = 0; j < size; j++)
		//		if (i >= 0 && j >= 0 && i < size && j < size)
		//		{
		//			if (i != 0) new_h[i, j] += (vh[i - 1, j] - vh[i, j]) * rate;
		//			if (i != size - 1) new_h[i, j] += (vh[i + 1, j] - vh[i, j]) * rate;
		//			if (j != 0) new_h[i, j] += (vh[i, j - 1] - vh[i, j]) * rate;
		//			if (j != size - 1) new_h[i, j] += (vh[i, j + 1] - vh[i, j]) * rate;
		//		}
		//Step 3
		//TODO: old_h <- h; h <- new_h;
		for (int i = 0; i < size; i++)
		{
			for (int j = 0; j < size; j++)
			{
				old_h[i, j] = h[i, j];
				h[i, j] = new_h[i, j];
			}
		}
		//Step 4: Water->Block coupling.
		//More TODO here.
	}
	

	// Update is called once per frame
	void Update () 
	{
		Mesh mesh = GetComponent ().mesh;
		Vector3[] X    = mesh.vertices;
		float[,] new_h = new float[size, size];
		float[,] h     = new float[size, size];


		//TODO: Load X.y into h.
		for (int i = 0; i < size; i++)
		{
			for (int j = 0; j < size; j++)
			{
				h[i, j] = X[i * size + j].y;
			}
		}
		if (Input.GetKeyDown ("r")) 
		{
			//TODO: Add random water.
			int i = Random.Range(0, size);
			int j = Random.Range(0, size);
			float rt = Random.value;
			float r = 0;
			for (int k = 0; k < 4; k++)
			{
				int ni = i + dirs[k, 0];
				int nj = j + dirs[k, 1];
				if (ni >= 0 && ni < size && nj >= 0 && nj < size)
				{
					r += h[ni, nj] * rt;
					h[ni, nj] -= h[ni, nj] * rt;
				}
			}
			//Debug.Log("add_r:" + r);
			h[i, j] += r;
		}
	
		for(int l=0; l<8; l++)
		{
			Shallow_Wave(old_h, h, new_h);
		}

		//TODO: Store h back into X.y and recalculate normal.
		for (int i = 0; i < size; i++)
		{
			for (int j = 0; j < size; j++)
			{
				X[i * size + j].y = h[i, j];
			}
		}

		mesh.vertices = X;
		mesh.RecalculateNormals();
	}
}

波动效果的公式

\[h_{i, j}^{\text {new }} \leftarrow h_{i, j}+\left(h_{i, j}-h_{i, j}^{\text {old }}\right) * \text { damping }+\left(h_{i-1, j}+h_{i+1, j}+h_{i, j-1}+h_{i, j+1}-4 h_{i, j}\right) * \text { rate. } \]

波动公式

经过转换其中的变量

这里简化ei=0 即 cube 块 直接和底面进行接触, 如果只考虑一个固体->流体的情况

vh 表示 虚拟高度

接下来就使用这张image

上面的图其实类似

只不过一个是真实高度一个是虚拟高度引起的江河湖海的高度变化.