fluid/solvers.edp at master · lvchien/fluid · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
// Solvers
macro NavierStokes()
{
    if (dc == 0)
    {
        solve NS ([ux, uy, p], [vx, vy, q], init = i) =
        int2d(Th)(
            alpha * (ux * vx + uy * vy)
            + nu * (dx(ux) * dx(vx) + dy(ux) * dy(vx) + dx(uy) * dx(vy) + dy(uy) * dy(vy))
            - p * q * 1e-8
            - p * dx(vx) - p * dy(vy)
            - dx(ux) * q - dy(uy) * q)
        - int2d(Th) (fx * vx + fy * vy)
        - int2d(Th) (alpha * convect([upx,upy], -dt, upx) * vx
                    + alpha * convect([upx,upy], -dt, upy) * vy)
        + on(1, ux = dcx1, uy = dcy1)
        + on(2, ux = dcx2, uy = dcy2)
        + on(3, ux = dcx3, uy = dcy3)
        + on(4, ux = dcx4, uy = dcy4)
        + on(20, ux = 0.0, uy = 0.0);
    }
    else if(dc == 1)
    {
        solve NS0 ([u0x, u0y, p0], [vx, vy, q], init = i) =
        int2d(Th)(
            alpha * (u0x * vx + u0y * vy)
            + (nu + sigma * Th.hmax) * (dx(u0x) * dx(vx) + dy(u0x) * dy(vx) + dx(u0y) * dx(vy) + dy(u0y) * dy(vy))
            - p0 * q * 1e-8
            - p0 * dx(vx) - p0 * dy(vy)
            - dx(u0x) * q - dy(u0y) * q
            )
        - int2d(Th) (fx * vx + fy * vy)
        - int2d(Th) (alpha * convect([upx,upy], -dt, up0x) * vx
                    + alpha * convect([upx,upy], -dt, up0y) * vy)
        + on(1, u0x = dcx1, u0y = dcy1)
        + on(2, u0x = dcx2, u0y = dcy2)
        + on(3, u0x = dcx3, u0y = dcy3)
        + on(4, u0x = dcx4, u0y = dcy4)
        + on(20, u0x = 0.0, u0y = 0.0);

    solve NS1 ([ux, uy, p], [vx, vy, q], init = i) =
        int2d(Th)(
            alpha * (ux * vx + uy * vy)
            + (nu + sigma * Th.hmax) * (dx(ux) * dx(vx) + dy(ux) * dy(vx) + dx(uy) * dx(vy) + dy(uy) * dy(vy))
            - p * q * 1e-8
            - p * dx(vx) - p * dy(vy)
            - dx(ux) * q - dy(uy) * q)
        - int2d(Th) (fx * vx + fy * vy)
        - int2d(Th) (alpha * convect([upx,upy], -dt, upx) * vx
                    + alpha * convect([upx,upy], -dt, upy) * vy)
        - int2d(Th) (
            sigma * Th.hmax * (dx(u0x) * dx(vx) + dy(u0x) * dy(vx) + dx(u0y) * dx(vy) + dy(u0y) * dy(vy))
            )
        + on(1, ux = dcx1, uy = dcy1)
        + on(2, ux = dcx2, uy = dcy2)
        + on(3, ux = dcx3, uy = dcy3)
        + on(4, ux = dcx4, uy = dcy4)
        + on(20, ux = 0.0, uy = 0.0);
    }
}//EOM

macro NaturalConvection()
{
    solve thermal (T, Tau, init = i) =
        int2d(Th) (alpha * T * Tau + dx(T) * dx(Tau) + dy(T) * dy(Tau))
        - int2d(Th) (Q * Tau)
        - int2d(Th) (alpha * convect([upx, upy], -dt, Tp) * Tau)
        + on(1, T = tbc1)
        + on(2, T = tbc2)
        + on(3, T = tbc3)
        + on(4, T = tbc4);
    fx = 0.0;
    fy = Ra * Pr * T;
    NavierStokes;
 } //EOM


macro Streamlines()
{
    solve strl(psi, phi) =
        int2d(Th)(dx(psi) * dx(phi) + dy(psi) * dy(phi))
        + int2d(Th)(-phi * (dy(ux) - dx(uy)))
        + on(1, 2, 3, 4, 10, 20, psi = 0);
} //EOM