xref: /libCEED/examples/fluids/postprocess/vortexshedding.py (revision ca567da4bd07bc8e9c3c34f09a40e154feb42909)
1import pandas
2import matplotlib.pyplot as plt
3import seaborn as sns
4import numpy as np
5import scipy.signal as sig
6
7
8def coeff(force, rho=1, u=1, D=1, zspan=0.2):
9    S = np.pi * D * zspan  # surface area
10    return 2 * force / (rho * u**2 * S)
11
12
13def shedding_period(df):
14    sample = df[df["Time"] > 70]  # once the initial transient has passed
15    peaks, _ = sig.find_peaks(sample["ForceY"])
16    period = np.diff(sample["Time"].iloc[peaks])
17    return period.mean()
18
19
20df = pandas.read_csv("force.csv")
21df["Drag Coefficient"] = coeff(df["ForceX"])
22df["Lift Coefficient"] = coeff(df["ForceY"])
23period = shedding_period(df)
24
25sns.set_theme(style="ticks")
26palette = sns.color_palette()
27fig, ax_drag = plt.subplots()
28ax_lift = ax_drag.twinx()
29
30sns.lineplot(data=df, x="Time", y="Drag Coefficient", ax=ax_drag, color=palette[0])
31sns.lineplot(data=df, x="Time", y="Lift Coefficient", ax=ax_lift, color=palette[1])
32ax_drag.set_title(f"Shedding period {period}")
33ax_drag.set_ylim(0.41, 0.49)
34ax_drag.tick_params(axis="y", colors=palette[0])
35ax_drag.yaxis.label.set_color(palette[0])
36ax_lift.tick_params(axis="y", colors=palette[1])
37ax_lift.yaxis.label.set_color(palette[1])
38
39# plt.savefig("vortexshedding.svg")
40plt.show()
41