numerical-programming/euler_2d.py at master · walterliu417/numerical-programming · 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
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
"""
Solves for potential flow in two dimensions by using continuity and the Navier-Stokes equations,
grossly simplified.
The flow is incompressible, inviscid (-> irrotational) and steady.
"""

"""
The governing equations are:
div(v) = 0                (continuity)
ρv • grad(v) = -grad(p)   (momentum balance (Navier-Stokes))
Since the flow is irrotational, we can write v = -grad(ϕ).
Then by continuity, -div(grad(ϕ)) = 0, i.e. Laplace's equation.
Finally, by Bernoulli's equation, p + 1/2 v^2 = const.
Hence the Bernoulli pressure = -(v • v) = -(grad(ϕ) • grad(ϕ)).
Imposing pressure boundary conditions will allow offseting to get the absolute pressure.
"""

import numpy as np
import time
import matplotlib.pyplot as plt
from copy import deepcopy as dcp
import random
from scipy import stats

import maths
from maths import *
from fem_base import *
from mesher import *

class PotentialFlowSolver2D:

    def __init__(self, bbox, wallfunc, V_in, poisson_func=lambda x, y: 0, resolution=5, tolerance=maths.tol):
        # The following BCs are used:
        # Inlet: -n • grad(ϕ) = V_in
        # Opening: ϕ = 0
        # Wall: -n • grad(ϕ) = 0
        # It is assumed that the left boundary is the inlet, other boundaries are openings,
        # and a no-slip wall is applied to the object with boundaries defined by wallfunc.

        # Generate initial mesh
        print("Starting meshing...\n")
        self.mesh2D = FluidMesh2D(bbox, wallfunc, resolution)
        self.V_in = V_in
        self.bbox = bbox
        self.tolerance = tolerance
        self.poisson_func = poisson_func


    def solve(self):
        # Optional to add the "driving function" f = -div(grad(ϕ)), for other applications.
        print("Solving...\n")
        start = time.time()
        # Form (n+1) x (n+1) stiffness matrix, A
        A = [[0 for j in range(self.mesh2D.num_nodes)] for i in range(self.mesh2D.num_nodes)]
        for element in self.mesh2D.elements:
            A[element.p0.id][element.p0.id] += element.area * (element.phi0.b ** 2 + element.phi0.c ** 2)
            A[element.p0.id][element.p1.id] += element.area * (element.phi0.b * element.phi1.b + element.phi0.c * element.phi1.c)
            A[element.p0.id][element.p2.id] += element.area * (element.phi0.b * element.phi2.b + element.phi0.c * element.phi2.c)
            A[element.p1.id][element.p0.id] += element.area * (element.phi1.b * element.phi0.b + element.phi1.c * element.phi0.c)
            A[element.p1.id][element.p1.id] += element.area * (element.phi1.b ** 2 + element.phi1.c ** 2)
            A[element.p1.id][element.p2.id] += element.area * (element.phi1.b * element.phi2.b + element.phi1.c * element.phi2.c)
            A[element.p2.id][element.p0.id] += element.area * (element.phi2.b * element.phi0.b + element.phi2.c * element.phi0.c)
            A[element.p2.id][element.p1.id] += element.area * (element.phi2.b * element.phi1.b + element.phi2.c * element.phi1.c)
            A[element.p2.id][element.p2.id] += element.area * (element.phi2.b ** 2 + element.phi2.c ** 2)
        A = np.array(A)

        # Form (n+1) x (n+1) residual matrix, R
        R = [[0 for j in range(self.mesh2D.num_nodes)] for i in range(self.mesh2D.num_nodes)]
        for edge in self.mesh2D.boundary_edges:
            if (edge.boundary == OPENING):
                R[edge.p1.id][edge.p1.id] += round((edge.length / 6) * (2/self.tolerance), -5)
                R[edge.p1.id][edge.p2.id] += round((edge.length / 6) * (1/self.tolerance), -5)
                R[edge.p2.id][edge.p1.id] += round((edge.length / 6) * (1/self.tolerance), -5)
                R[edge.p2.id][edge.p2.id] += round((edge.length / 6) * (2/self.tolerance), -5)
        R = np.array(R)

        # Form (n+1) load vector, b
        b = [0 for i in range(self.mesh2D.num_nodes)]
        for element in self.mesh2D.elements:
            b[element.p0.id] += (element.area / 3) * self.poisson_func(element.p0.x, element.p0.y)
            b[element.p1.id] += (element.area / 3) * self.poisson_func(element.p1.x, element.p1.y)
            b[element.p2.id] += (element.area / 3) * self.poisson_func(element.p2.x, element.p2.y)
        b = np.array(b)

        # Form (n+1) residual vector, r
        # Note vector b = int(f * v dK) = 0, since f = 0 (special case of Poissons')
        r = [0 for i in range(self.mesh2D.num_nodes)]
        for edge in self.mesh2D.boundary_edges:
            if edge.boundary == INLET:
                r[edge.p1.id] += (edge.length / 2) * self.V_in
                r[edge.p2.id] += (edge.length / 2) * self.V_in
            elif edge.boundary == OUTLET:
                r[edge.p1.id] += (edge.length / 2) * -self.V_in
                r[edge.p2.id] += (edge.length / 2) * -self.V_in
        r = np.array(r)

        # Symmetry check
        try:
            assert np.allclose((A+R), (A+R).T)
        except:
            raise Warning("Stiffness matrix is not symmetric??")

        # Solve the linear system
        linstart = time.time()
        self.coeffs = np.linalg.solve(A+R, b+r)
        print(f"Linear solution completed in {time.time() - start} seconds.")
        self.solution = FEMSolution2D(self.mesh2D.points, self.mesh2D.elements, self.coeffs)

        x, y, z = [], [], []
        real = []
        for p in self.mesh2D.points:
            x.append(p.x)
            y.append(p.y)
            #z.append(self.solution(p.x, p.y))
            r = float((p.x * p.x + p.y * p.y) ** 0.5)
            try:
                if p.x > 0:
                    theta = np.arctan(float(p.y) / float(p.x))
                elif p.x < 0 and p.y >= 0:
                    theta = np.arctan(float(p.y) / float(p.x)) + np.pi
                elif p.x < 0 and p.y < 0:
                    theta = np.arctan(float(p.y) / float(p.x)) - np.pi

            except ZeroDivisionError:
                if p.y > 0: theta = np.pi / 2
                else: p.y = -np.pi / 2
            real.append(-self.V_in * r * (1 + 16 / (r ** 2)) * np.cos(theta))
        print("Solution completed in", time.time() - start, "seconds.")
        print("Velocity: ", self.V_in)
        fig, (ax1, ax2) = plt.subplots(2)
        c1 = ax1.tricontourf(x, y, [[e.p0.id, e.p1.id, e.p2.id] for e in self.mesh2D.elements], self.coeffs, 40)
        c2 = ax2.tricontourf(x, y, [[e.p0.id, e.p1.id, e.p2.id] for e in self.mesh2D.elements], real, 40)
        plt.colorbar(c1, ax=ax1)
        plt.colorbar(c2, ax=ax2)
        plt.show()
        return self.solution


# Testing area
if __name__ == "__main__":
    def bus(x, y):
        if (-4 < x < 4) and (-2 < y < 0.5 + np.sin(float(4 * x))):
            return False
        return True

    def cylinder(x, y):
        return x ** 2 + y ** 2 >= 16

    def source(x, y):
        if x < -10 + tol:
            return 1
        return 0

    flow_speed = 10

    solution = PotentialFlowSolver2D((-10, -10, 10, 10), cylinder, flow_speed, resolution=70)
    result = solution.solve()
    result.plot_velocity_and_pressure(flow_speed)