2.txt

圣维南方程组的显示差分格式源代码

using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Text;
using System.Windows.Forms;
using System.IO;

namespace 显式差分
{
    public partial class Form1 : Form
    {
        public Form1()
        {
            InitializeComponent();
        }

        private void button1_Click(object sender, EventArgs e)
        {
            double[,] zz = new double[500, 11]; double[,] qq = new double[500, 11];
            double[] z = new double[11] { 15.518,15.515,15.51,15.51,15.51,15.51,15.51,15.5,15.5,15.5,15.5};
            double[] q = new double[11]; double[] qp = new double[11]; double[] zp = new double[11];
            double[] z0 = new double[11] { 11.018,10.915,10.81,10.71,10.61,10.51,10.41,10.3,10.2,10.1,10.0};
            double[] ds = new double[11] { 513.4,512.1,510.9,509.9,509.0,508.1,507.4,506.7,506.1,505.6,0.0};
            double dt = 60; int nf = 11; double g = 9.8; double cn = 0.013; double s = 0.0002; double cm = 3; int j = 0; double b = 5;
            int i;
            double h;
            double bs;
            double a;
            double p;
            double r;
            double w0;
            double w1;
            double cn0;
            double zl=0,ql=0,zr,qr;
            FileStream fs = new FileStream(Application.StartupPath + @"\output.xls", FileMode.OpenOrCreate);
            StreamWriter sw = new StreamWriter(fs, Encoding.Unicode);

            for (i = 0; i < nf; i++)
            {
                q[i] = 30;
                qq[0, i] = 30;
                zz[0, i] = z[i];
            }
            int t = 0;
            
            while (t <= 360)
            {

                j = j + 1;
                t = (int)(dt / 60 * j);
                qp[nf - 1] = 30 + 120 / 20 * t; zp[0] = 15.518;
                if(qp[nf-1]>150)
                    qp[nf-1]=150;
                for (i = 0; i < nf; i++)
                {
                    h = z[i] - z0[i];
                    bs = b + 2 * cm * h;
                    a = (b + cm * h) * h;
                    p = b + Math.Sqrt(1 + cm * cm) * h * 2;
                    r = a / p;
                    w0 = bs * (q[i] / a - Math.Sqrt(g * a / bs));
                    w1 = bs * (q[i] / a + Math.Sqrt(g * a / bs));
                    cn0 = -bs * s * Math.Pow((q[i] / a), 2) + g * q[i] * q[i] * cn * cn / a / Math.Pow(r, 4.0 / 3.0);
                    if (i != 0)
                    {
                        zl = z[i] + dt / ds[i - 1] * (z[i - 1] - z[i]) * w1 / bs;
                        ql = q[i] + dt / ds[i - 1] * (q[i - 1] - q[i]) * w1 / bs;
                        if (i == nf - 1)
                        {
                            zp[nf - 1] = zl + (cn0 * dt + qp[nf - 1] - ql) / w0;
                        }
                    }
                    if (i != nf - 1)
                    {
                        zr = z[i] + dt / ds[i] * (z[i] - z[i + 1]) * w0 / bs;
                        qr = q[i] + dt / ds[i] * (q[i] - q[i + 1]) * w0 / bs;
                        if (i != 0)
                            zp[i] = (qr - ql + w0 * zl - w1 * zr) / (w0 - w1);
                        qp[i] = qr + w1 * (zp[i] - zr) - cn0 * dt;
                    }
                    sw.WriteLine(t.ToString() + "\t" + i.ToString() + "\t" + qp[i].ToString() + "\t" + zp[i].ToString());

                }
                for (i = 0; i < nf; i++)
                {
                    q[i] = qp[i];
                    qq[j, i] = qp[i];
                    z[i] = zp[i];
                    zz[j, i] = zp[i];
                }
                sw.WriteLine(t.ToString() + "\t" + q[0].ToString() + "\t" + z[0].ToString() + "\t" + q[nf - 1].ToString() + "\t" + z[nf - 1].ToString());
            }

        }
    }
}