Bibliographie :
Idées :
import pandas as pd
from io import StringIO
import holoviews as hv
hv.extension('bokeh')
import numpy as np
Les équations du potentiel :
$$ \vec{v} = \vec{grad} \, \varphi $$$$ \Delta \varphi = 0 $$Conditions aux limites :
En amont : $$ \vec{v} = \vec{v_0} $$
Et autour de l'obstacle : $$ \vec{v} \cdot \vec{n} = 0 $$
Il faut aussi réfléchir aux conditions sur les bords du maillage : est-ce qu'on force la vitesse à être nulle ? (hypothèse de zone non perturbée par l'écoulement ?)
On définit un motif d'aile avec des points discrets (ce sera plus simple à interpoler plus tard).
mesh = """315.3,242.75
317.2142857142857,234.85714285714286
317.2142857142857,249.85714285714286
320.5,254.83333333333334
322.22727272727275,230.54545454545453
325.5,257.625
328,227.25
328.5,261.5625
334,224.5
333.5,263.3888888888889
337.5,267.05555555555554
340.25,222.25
342.5,268.5
345.92857142857144,220.78571428571428
346.2142857142857,271.57142857142856
352.2142857142857,220.07142857142858
351.625,273
356.4,276.3
358,219
362.05555555555554,277.94444444444446
364,219
367.3,280.65
370,219
372.72222222222223,282.1666666666667
376,219
378.27777777777777,284.77777777777777
381.7857142857143,220.07142857142858
383.625,286
388,220
389.27777777777777,288.77777777777777
394,220.66666666666666
394.72222222222223,290.1666666666667
400,221.5
400.05555555555554,292.3333333333333
406,222.58333333333334
406,293.625
412,223.25
412,295.625
418,224.33333333333334
418.2,297.9
424,225.58333333333334
424.25,299.375
430.07142857142856,226.78571428571428
430.07142857142856,300.7857142857143
435.92857142857144,227.92857142857142
435.875,302.75
441.75,229.875
442.125,304.5
448,231.625
448.07142857142856,305.7857142857143
454,233.625
453.875,307.75
460.2,235.9
459.92857142857144,308.92857142857144
466,237.625
465.75,310.875
472,239.625
472.125,312.5
478,241.625
478.07142857142856,313.7857142857143
483.8,244.1
484.125,315.5
489.75,245.875
490.2142857142857,316.64285714285717
495.1666666666667,248.55555555555554
495.92857142857144,317.92857142857144
500.72222222222223,250.16666666666666
501.875,319.75
506.27777777777777,252.77777777777777
508.07142857142856,320.7857142857143
511.94444444444446,254.72222222222223
514.125,322.5
517.6111111111111,256.94444444444446
520.2142857142857,323.64285714285717
523.2777777777778,259.77777777777777
525.9285714285714,324.92857142857144
528.9444444444445,261.72222222222223
532.125,326.5
534.6,264.2
538.2142857142857,327.64285714285717
540.3888888888889,267
544.0714285714286,328.7857142857143
546.0555555555555,268.94444444444446
549.7857142857143,330.07142857142856
551.1666666666666,271.55555555555554
556,331.625
557,274
562,332.3333333333333
563.125,276.5
568,333.5
569.375,279.25
574,334.6666666666667
575,282
580.1666666666666,284.55555555555554
580,336.625
586,287
586.2142857142857,337.64285714285717
591.1666666666666,289.55555555555554
592.2142857142857,338.64285714285717
596.5,291.3888888888889
598.0714285714286,339.7857142857143
600.5,294.55555555555554
603.9285714285714,340.92857142857144
605.875,295.75
609.7857142857143,342.07142857142856
611.1666666666666,298.8888888888889
615.6428571428571,343.2142857142857
616.6,301.2
621.4,304.3
622,344.625
626.6,306.2
628.3571428571429,345.5
632.125,308.8125
634.2142857142857,346.64285714285717
637.625,311
640.2142857142857,347.64285714285717
642.4,314.3
646.0714285714286,348.7857142857143
647.6,316.2
652.4,319.3
652.0714285714286,349.7857142857143
657.6,321.2
657.9285714285714,350.92857142857144
662.4,324.3
663.9285714285714,351.92857142857144
667.6,326.2
669.9285714285714,352.92857142857144
672.4,329.3
676,353.6666666666667
677.6,331.2
682.4,334.3
682,354.6666666666667
687.6428571428571,336.5
688,355.6666666666667
694,339.75
693.7857142857143,357.07142857142856
700,342.75
699.7857142857143,358.07142857142856
706,345.9166666666667
706,358.6666666666667
712,349.3333333333333
711.7857142857143,360.07142857142856
718,351.75
717.7857142857143,361.07142857142856
724,354.1666666666667
724.4,362.25
730,356.5
729.625,363.0625
736,359.8333333333333
736,364.6666666666667
741.0714285714286,364.07142857142856
"""
df = pd.read_csv(StringIO(mesh), header=None)
df.columns = ['x', 'y']
df = df.sort_values(by='x')
df.head()
x | y | |
---|---|---|
0 | 315.300000 | 242.750000 |
1 | 317.214286 | 234.857143 |
2 | 317.214286 | 249.857143 |
3 | 320.500000 | 254.833333 |
4 | 322.227273 | 230.545455 |
hv.Scatter(df, 'x', 'y')
On va trier les points par plus proche voisin pour avoir un contour orienté correctement.
from scipy.spatial import KDTree
def reorder_mesh(xy_points):
"""Reorders mesh points by nearest neighbor."""
tree = KDTree(xy_points)
next_point = 0
ordered_mesh = [next_point]
while len(ordered_mesh) < xy_points.shape[0]:
dists, indices = tree.query(xy_points[ordered_mesh[-1]], k=20)
next_point = [index for index in indices if index not in ordered_mesh]
if len(next_point) > 0:
ordered_mesh.append(next_point[0])
else:
break
ordered_xy = xy_points[ordered_mesh]
return ordered_xy
df = pd.DataFrame(reorder_mesh(df.values), columns=['x', 'y'])
poly_mesh = hv.Polygons([{('x', 'y'): df.values, 'level':0}])
poly_mesh
On va choisir une grille qui englobe l'aile définie précédemment.
dx = dy = 10
x = np.arange(0, 1000, dx)
y = np.arange(0, 400, dy)
X, Y = np.meshgrid(x, y)
bounds = (x.min(), y.min(), x.max(), y.max())
hv.Image(X, bounds=bounds, vdims='xval') + hv.Image(Y, bounds=bounds, vdims='yval')
On va se ramener à la résolution du système linéaire $A \varphi = b$. Pour cela on distinguer d'abord les points du maillage sur lesquels on va résoudre le potentiel et qui ne sont pas à l'intérieur du profil.
On va construire une grille booléenne pour ce faire :
from shapely.geometry import Point
from shapely.geometry.polygon import Polygon
bool_mask = []
polygon = Polygon(df.values)
for xy in zip(X.flatten(), Y.flatten()):
point = Point(xy)
bool_mask.append(polygon.contains(point))
bool_mask = np.array(bool_mask).reshape(X.shape)
hv.Image(bool_mask, bounds=bounds)
np.sum(bool_mask)
212
On peut maintenant boucler sur tous les points du maillage dans l'ordre qui nous intéresse. Il va nous falloir à la fois un mapping rapide entre les coordonnées d'un point et son indice dans le tableau des inconnues.
mesh_ids = np.nonzero(~bool_mask)
id_to_coords = {ind:(r, c) for ind, (r, c) in enumerate(zip(mesh_ids[0], mesh_ids[1]))}
coords_to_id = dict(zip(id_to_coords.values(), id_to_coords.keys()))
On crée une grande matrice et on la remplit :
N = len(id_to_coords)
A = np.zeros((N, N))
A.shape
R, C = X.shape
alpha = 1/(2. * dx)
# "normal" cell
for nid in id_to_coords:
r, c = id_to_coords[nid]
if all((rr, cc) in coords_to_id for rr, cc in [(r, c),
(r+1, c),
(r-1, c),
(r, c+1),
(r, c-1)]):
if r>0 and r<R-1 and c>0 and c<C-1:
A[nid, nid] += -4 * alpha
A[nid, coords_to_id[r+1, c]] += 1 * alpha
A[nid, coords_to_id[r-1, c]] += 1 * alpha
A[nid, coords_to_id[r, c+1]] += 1 * alpha
A[nid, coords_to_id[r, c-1]] += 1 * alpha
On peut maintenant rajouter les conditions sur le bord où on fixe la vitesse (conditions de Neumann donc vu qu'on est en potentiel).
b = np.zeros((N,))
v0 = 1.
c = 0
for r in range(R):
nid = coords_to_id[(r, c)]
r, c = id_to_coords[nid]
if all((rr, cc) in coords_to_id for rr, cc in [(r, c),
(r, c+1)]):
A[nid, coords_to_id[r, c+1]] += 1
A[nid, coords_to_id[r, c]] -= 1
b[nid] = v0 * dx
Finalement, on peut rajouter les conditions de non-pénétration sur le bord de l'obstacle.
Et maintenant on peut résoudre :
%%opts Scatter [width=800]
hv.Scatter(np.sum(np.abs(A), axis=1))
On peut dessiner les lignes qui sont nulles par leur localisation spatiale.
zero_indices = np.where(np.sum(np.abs(A), axis=1)==0)[0]
id_to_coords[0]
(0, 0)
flat_XY = np.c_[Y.flatten(), X.flatten()]
zero_pts = np.array([flat_XY[np.array(id_to_coords[nid])] for nid in zero_indices])
hv.Scatter(zero_pts)
hv.Image(bool_mask, bounds=bounds).options(alpha=.5) * hv.Scatter(zero_pts)
np.linalg.solve(A, b)
--------------------------------------------------------------------------- LinAlgError Traceback (most recent call last) <ipython-input-28-eb79222c8972> in <module>() ----> 1 np.linalg.solve(A, b) C:\Anaconda3\lib\site-packages\numpy\linalg\linalg.py in solve(a, b) 382 signature = 'DD->D' if isComplexType(t) else 'dd->d' 383 extobj = get_linalg_error_extobj(_raise_linalgerror_singular) --> 384 r = gufunc(a, b, signature=signature, extobj=extobj) 385 386 return wrap(r.astype(result_t, copy=False)) C:\Anaconda3\lib\site-packages\numpy\linalg\linalg.py in _raise_linalgerror_singular(err, flag) 88 89 def _raise_linalgerror_singular(err, flag): ---> 90 raise LinAlgError("Singular matrix") 91 92 def _raise_linalgerror_nonposdef(err, flag): LinAlgError: Singular matrix
Si la matrice est trop singulière, on peut toujours utiliser la pseudo inverse...
np.linalg.lstsq(A, b)