Skip to content

Commit 5f991e3

Browse files
Merge pull request #148 from KratosMultiphysics/fsi_mok_iga_fem
[CoSimulationApplication] FSI-Mok Benchmark IGA-FEM
2 parents 85b3079 + 10cb4ef commit 5f991e3

16 files changed

+21375
-0
lines changed

co_simulation/validation/README.md

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,8 @@ This folder contains the validation cases:
44

55
- [Mok FSI benchmark](https://github.com/KratosMultiphysics/Examples/blob/master/co_simulation/validation/fsi_mok/README.md)
66

7+
- [Mok FSI benchmark IGA-FEM](https://github.com/KratosMultiphysics/Examples/blob/master/co_simulation/validation/fsi_mok_iga_fem/README.md)
8+
79
- [Dam break against a flexible wall FSI benchmark](https://github.com/KratosMultiphysics/Examples/blob/master/co_simulation/validation/dam_break_flex_wall/README.md)
810

911
- [FSI Turek Benchmark](https://github.com/KratosMultiphysics/Examples/tree/master/co_simulation/validation/fsi_turek_FSI2)
Lines changed: 66 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,66 @@
1+
# Mok FSI benchmark
2+
3+
**Author:** [Juan Ignacio Camarotti](https://github.com/juancamarotti)
4+
5+
**Kratos version:** 10.2
6+
7+
**Source files:** [FSI-Mok IGA-FEM](https://github.com/KratosMultiphysics/Examples/tree/master/co_simulation/validation/fsi_mok_iga_fem/source)
8+
9+
## Case Specification
10+
11+
This is a 2D FSI simulation of the Mok benchmark test. It consists in a 2D convergent fluid channel that contains a flexible wall structure attached to its bottom wall. The main challenge of the test is that the densities of the fluid and the structure have a similar order of magnitude, leading to a strongly coupled problem in where large interaction between the two fields appears. The reference solutions have been taken from Mok (2001) and Valdés (2007). The following applications of Kratos are used:
12+
* CoSimulationApplication
13+
* MappingApplication
14+
* MeshMovingApplication
15+
* FluidDynamicsApplication
16+
* StructuralMechanicsApplication
17+
* IgaApplication
18+
* LinearSolversApplication
19+
20+
The problem geometry as well as the boundary conditions are sketched below.
21+
<p align="center">
22+
<img src="data/Mok_benchmark_geometry.png" alt="Mok benchmark geometry." style="width: 600px;"/>
23+
</p>
24+
25+
Regarding the inlet velocity, the next parabolic profile is imposed
26+
27+
<p align="center">
28+
<img src="data/Mok_inlet_formula_1.png" alt="Mok inlet profile." style="width: 200px;"/>
29+
</p>
30+
31+
where the time dependent reference velocity is defined as
32+
33+
<p align="center">
34+
<img src="data/Mok_inlet_formula_2.png" alt="Mok velocity formula." style="width: 200px;"/>
35+
</p>
36+
37+
A Newtonian constitutive law is considered in the fluid domain. The fluid characteristic parameters are:
38+
* Density (&rho;): 956 _Kg/m<sup>3</sup>_
39+
* Kinematic viscosity (&nu;): 0.145 _m<sup>2</sup>/s_
40+
41+
On the other hand, a linear elastic plane stress constitutive law with unit thickness is considered in the structure domain. The structure characteristic parameters are
42+
* Density (&rho;): 1500 _Kg/m<sup>3</sup>_
43+
* Elastic modulus (E): 2.30000E+06 _Pa_
44+
* Poisson ratio (&nu;): 0.45
45+
46+
The time step is 0.1 seconds, while the total simulation time is 25.0 seconds.
47+
48+
The mesh was created with the [KratosSalomePlugin](https://github.com/KratosMultiphysics/KratosSalomePlugin/tree/master/tui_examples/mok_fsi). Check this example which can be easily adapted to different mesh sizes.
49+
50+
The mapping strategy employed in this simulation is the nearest neighbor mapper, specifically tailored for partitioned IGA-FEM simulations.
51+
52+
## Results
53+
The structural domain for the problem described above was discretized using an IBRA mesh with Kirchhoff-Love shell elements (3-parameter formulation). For the fluid domain, a mesh consisting of approximately 6000 linear triangular elements was employed. The resulting velocity field, along with the deformed geometry, is presented below.
54+
55+
<p align="center">
56+
<img src="data/flow_field_vel_t25.png" alt="Obtained velocity field (t = 25.0)." style="width: 600px;"/>
57+
</p>
58+
59+
<p align="center">
60+
<img src="data/Mok_ux.png" alt="Point A horizontal displacement comparison." style="width: 600px;"/>
61+
</p>
62+
63+
## References
64+
D.P. Mok. Partitionierte Lösungsansätze in der Strukturdynamik und der Fluid−Struktur−Interaktion. PhD thesis: Institut für Baustatik, Universität Stuttgart, 2001. [http://dx.doi.org/10.18419/opus-147](http://dx.doi.org/10.18419/opus-147)
65+
66+
G. Valdés. Nonlinear Analysis of Orthotropic Membrane and Shell Structures Including Fluid-Structure Interaction. PhD thesis: Universitat Politècnica de Catalunya, 2007. [http://www.tdx.cat/handle/10803/6866](http://www.tdx.cat/handle/10803/6866)
157 KB
Loading
2.18 KB
Loading
8.17 KB
Loading
75.4 KB
Loading
89.8 KB
Loading
Lines changed: 117 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,117 @@
1+
import KratosMultiphysics as KM
2+
import KratosMultiphysics.IgaApplication
3+
from KratosMultiphysics.CoSimulationApplication.co_simulation_analysis import CoSimulationAnalysis
4+
5+
# External imports
6+
import importlib
7+
import matplotlib.pyplot as plt
8+
import numpy as np
9+
from scipy.interpolate import griddata
10+
11+
import sys
12+
import time
13+
14+
class CustomCoSimulationAnalysis(CoSimulationAnalysis):
15+
def __init__(self, parameters):
16+
"""Call the base class constructor."""
17+
super().__init__(parameters)
18+
self.time_history = []
19+
self.disp_x_history_A = []
20+
self.disp_x_history_B = []
21+
22+
def OutputSolutionStep(self):
23+
super().OutputSolutionStep()
24+
25+
solver = self._GetSolver()
26+
27+
if hasattr(solver, "model") and hasattr(solver.model, "solver_wrappers"):
28+
structure_solver = solver.model.solver_wrappers.get("structure")
29+
if structure_solver and hasattr(structure_solver, "model"):
30+
sub_model = structure_solver.model
31+
model_part_name = "IgaModelPart"
32+
33+
if hasattr(sub_model, "HasModelPart") and sub_model.HasModelPart(model_part_name):
34+
model_part = sub_model.GetModelPart(model_part_name)
35+
36+
wet_interface_sub_model_part = model_part.GetSubModelPart("Load_4")
37+
current_time = wet_interface_sub_model_part.ProcessInfo[KratosMultiphysics.TIME]
38+
39+
# Only do this the first time
40+
if not hasattr(self, "tracking_conditions_initialized"):
41+
self.tracking_conditions_initialized = True
42+
43+
# Choose target locations
44+
self.point_A = (0.50025, 0.25) # Replace with your target coordinates
45+
self.point_B = (0.50, 0.125)
46+
47+
def get_closest_condition(target_point):
48+
min_dist = float("inf")
49+
closest_cond = None
50+
for cond in wet_interface_sub_model_part.Conditions:
51+
center = cond.GetGeometry().Center()
52+
dist = (center.X - target_point[0])**2 + (center.Y - target_point[1])**2
53+
if dist < min_dist:
54+
min_dist = dist
55+
closest_cond = cond
56+
return closest_cond
57+
58+
self.condition_A = get_closest_condition(self.point_A)
59+
self.condition_B = get_closest_condition(self.point_B)
60+
61+
# --- Evaluate displacement at A ---
62+
condition = self.condition_A
63+
geom = condition.GetGeometry()
64+
N = geom.ShapeFunctionsValues()
65+
66+
solution_A = 0.0
67+
for i, node in enumerate(condition.GetNodes()):
68+
solution_A += node.GetSolutionStepValue(KratosMultiphysics.DISPLACEMENT_X) * N[0, i]
69+
70+
self.time_history.append(current_time)
71+
self.disp_x_history_A.append(solution_A)
72+
73+
# --- Evaluate displacement at B ---
74+
condition = self.condition_B
75+
geom = condition.GetGeometry()
76+
N = geom.ShapeFunctionsValues()
77+
78+
solution_B = 0.0
79+
for i, node in enumerate(condition.GetNodes()):
80+
solution_B += node.GetSolutionStepValue(KratosMultiphysics.DISPLACEMENT_X) * N[0, i]
81+
82+
self.disp_x_history_B.append(solution_B)
83+
84+
def Finalize(self):
85+
super().Finalize()
86+
87+
import matplotlib.pyplot as plt
88+
89+
# Enable LaTeX-style rendering
90+
plt.rcParams["text.usetex"] = True
91+
plt.rcParams["font.family"] = "serif"
92+
93+
# Plot
94+
plt.figure()
95+
plt.plot(self.time_history, self.disp_x_history_A, label=r"$u_x$ at Point A")
96+
plt.plot(self.time_history, self.disp_x_history_B, label=r"$u_x$ at Point B")
97+
98+
# Labels with LaTeX formatting
99+
plt.xlabel(r"\textbf{Time} [s]")
100+
plt.ylabel(r"\textbf{Displacement} $u_x$ [m]")
101+
plt.title(r"\textbf{Displacement vs Time (Nearest Neighbor Mapper)}", fontsize=14)
102+
103+
# Grid and legend
104+
plt.grid(True)
105+
plt.legend(loc="best")
106+
107+
# Save and show
108+
plt.savefig("displacement_vs_time.pdf", dpi=300, bbox_inches="tight")
109+
plt.show()
110+
111+
if __name__ == "__main__":
112+
113+
with open("ProjectParametersCoSim.json", 'r') as parameter_file:
114+
parameters = KratosMultiphysics.Parameters(parameter_file.read())
115+
116+
simulation = CustomCoSimulationAnalysis(parameters)
117+
simulation.Run()

0 commit comments

Comments
 (0)