Dem rotating drum example post processing (#1394)

Description
Added the post processing for the rotating drum example.

CHANGELOG.md

@@ -3,6 +3,11 @@
 All notable changes to the Lethe project will be documented in this file.
 The format is based on [Keep a Changelog](http://keepachangelog.com/).
 
+## [Master] - 2024-12-16
+
+### Added
+
+- MINOR A new post-processing code is now available for the rotating drum example. This code output the average velocity profile of particles perpendicular to the free surface. 
 
 ## [Master] - 2024-12-03
 

examples/dem/3d-rotating-drum/data_exp.csv

@@ -0,0 +1,36 @@
+-0.3480331262939959, -0.03030549898167012
+-0.33105590062111806, -0.03494908350305502
+-0.31656314699792965, -0.03262729124236252
+-0.29668737060041406, -0.03714867617107946
+-0.2788819875776397, -0.028228105906313666
+-0.2544513457556936, -0.0413034623217923
+-0.22339544513457552, -0.03898167006109983
+-0.20890269151138713, -0.04350305498981673
+-0.19523809523809524, -0.04496945010183302
+-0.1546583850931677, -0.04949083503054992
+-0.143064182194617, -0.04802443991853364
+-0.12318840579710147, -0.05169042769857439
+-0.09089026915113874, -0.05389002036659882
+-0.04120082815734988, -0.05572301425661917
+-0.01345755693581785, -0.05865580448065178
+0.013871635610766042, -0.06012219959266805
+0.015942028985507284, -0.06219959266802446
+0.039544513457556996, -0.06513238289205708
+0.05486542443064191, -0.06708757637474547
+0.07101449275362326, -0.06953156822810595
+0.08260869565217382, -0.0709979633401222
+0.09171842650103518, -0.07405295315682286
+0.0966873706004141, -0.07551934826883912
+0.1049689440993789, -0.07747454175152754
+0.10828157349896483, -0.07967413441955197
+0.11821946169772257, -0.0818737270875764
+0.11987577639751551, -0.08431771894093691
+0.13022774327122144, -0.08566191446028515
+0.12194616977225692, -0.08883910386965375
+0.12815734989648037, -0.09079429735234221
+0.13105590062111805, -0.09543788187372712
+0.13561076604554861, -0.09861507128309574
+0.13602484472049686, -0.10118126272912428
+0.13685300207039341, -0.10399185336048881
+0.1372670807453416, -0.10680244399185342
+0.14140786749482398, -0.11303462321792265

examples/dem/3d-rotating-drum/post_processing_rotating_drum.py

@@ -0,0 +1,109 @@
+# SPDX-FileCopyrightText: Copyright (c) 2024 The Lethe Authors
+# SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
+
+import numpy as np
+import matplotlib.pyplot as plt
+import pandas as pd
+import argparse
+from lethe_pyvista_tools import *
+from numpy.ma.core import arctan
+from numpy.ma.extras import average
+
+parser = argparse.ArgumentParser(description='Arguments for calculation of the velocity profile in the rotating drum')
+parser.add_argument("-f", "--folder", type=str, help="Folder path", required=True)
+parser.add_argument("-i", "--input", type=str, help="Name of the input file", required=True)
+args, leftovers=parser.parse_known_args()
+
+folder=args.folder
+filename=args.input
+
+# Experimental date extracted using Web plot digitizer
+df_x_zero_exp = pd.read_csv(folder + "/" +filename, header=None,
+                            names=['Column1', 'Column2'])
+
+# Starting vtu id. Here, we are interested in the last 100 vtu
+start = 900
+
+# Load lethe data
+pvd_name = 'out.pvd'
+ignore_data = ['type', 'volumetric contribution', 'torque', 'fem_torque',
+               'fem_force']
+particle = lethe_pyvista_tools("./", "rotating-drum.prm",
+                               pvd_name, ignore_data=ignore_data)
+time = particle.time_list
+
+# We fix the angle at 27 deg like is was mentioned in the experimental article
+angle = 27. * 2. * np.pi / 360.
+
+# H is the normal distance between the free surface and the center of the cylinder
+H = -0.03
+
+# Point where the local averages are made
+y_graph = np.linspace(-0.12, H, 100) # -0.12 is the cylinder radius
+vel_value_x = np.zeros_like(y_graph, float)
+vel_value_y = np.zeros_like(y_graph, float)
+x_graph = 0.
+D = 1. * 0.003 # D is the distance in which particle are considered for the local average
+
+# Change frame of reference. (loc is for local)
+# The X direction is parallel to the free surface and pointing down.
+x_loc_unit = np.array([-np.cos(angle), -np.sin(angle)])
+y_loc_unit = np.array([-np.sin(angle), np.cos(angle)])
+
+for i in range(start, len(time)):
+    print("Time step: ", i)
+    df_load = particle.get_df(i)
+
+    # We do a scalar product to find the velocity in the new frame of reference.
+    df = pd.DataFrame(df_load.points, columns=['x', 'y', 'z'])
+    df['v_x_loc'] = pd.DataFrame(df_load['velocity'][:, 1]) * x_loc_unit[
+        0] + pd.DataFrame(df_load['velocity'][:, 2]) * x_loc_unit[1]
+
+    # Same here, scalar product
+    df['v_y_loc'] = pd.DataFrame(df_load['velocity'][:, 1]) * y_loc_unit[
+        0] + pd.DataFrame(df_load['velocity'][:, 2]) * y_loc_unit[1]
+
+
+    # The cylinder is 0.12 m long.
+    # To reduce the wall effect in our calculation, we only focus on the particles
+    # 0.02 m from the front and back walls
+    cut_off = 0.10
+    condition = (df['x'] < cut_off) & (df['x'] > -cut_off)
+    df = df[condition]
+
+    df['x_loc'] = (df['y'] ) * x_loc_unit[0] + (df['z']) * \
+                  x_loc_unit[1]
+    df['y_loc'] = (df['y'] ) * y_loc_unit[0] + (df['z']) * \
+                  y_loc_unit[1]
+
+    df_filtered = df.copy()
+
+    for index, j in enumerate(y_graph):
+        # Compute the distance between each particle and the local average calculation point.
+        df_filtered['dist'] = (df_filtered['x_loc'] ** 2. + (df_filtered['y_loc'] - j) ** 2.) ** 0.5
+        condition_2 = (df_filtered['dist'] < D)
+        df_filtered_filtered = df_filtered[condition_2]
+
+        # Number of particle in the local average zone
+        n = len(df_filtered_filtered)
+
+        # Average velocity in the local average
+        tot_vel_x = df_filtered_filtered['v_x_loc'].sum()
+        vel_value_x[index] += tot_vel_x / n
+        tot_vel_y = df_filtered_filtered['v_y_loc'].sum()
+        vel_value_y[index] += tot_vel_y / n
+
+# Average over time of the averages
+vel_value_x = vel_value_x / (len(time) - start)
+vel_value_y = vel_value_y / (len(time) - start)
+
+plt.plot(vel_value_x, y_graph, label="Lethe")
+plt.plot(-df_x_zero_exp["Column1"].to_numpy(),
+         df_x_zero_exp["Column2"].to_numpy(), "xk", label="RPT")
+plt.xlabel("u (m/s)")
+plt.ylabel("y (m)")
+plt.plot(y_graph * 1.214749159, y_graph, "-.", label=r"$\omega r$")
+plt.legend()
+plt.title(f"Average velocity parallel to the free surface \ndepending on the depth from the center of the cylinder")
+plt.show()
+

examples/dem/3d-rotating-drum/rotating-drum.prm

@@ -34,6 +34,7 @@ subsection model parameters
     set frequency           = 20000
   end
   set particle particle contact force method = hertz_mindlin_limit_overlap
+  set rolling resistance torque method       = constant_resistance
   set particle wall contact force method     = nonlinear
   set integration method                     = velocity_verlet
 end
@@ -53,12 +54,14 @@ subsection lagrangian physical properties
     set young modulus particles           = 1e7
     set poisson ratio particles           = 0.24
     set restitution coefficient particles = 0.97
-    set friction coefficient particles    = 0.3
+    set friction coefficient particles    = 0.85
+    set rolling friction particles        = 0.025
   end
   set young modulus wall           = 1e7
   set poisson ratio wall           = 0.24
   set restitution coefficient wall = 0.85
-  set friction coefficient wall    = 0.35
+  set friction coefficient wall    = 0.85
+  set rolling friction wall        = 0.025
 end
 
 subsection restart