<jupyter_start><jupyter_text>Flow in a nano-porous materialIn this notebook we'll explore flow induced in a nano-porous material. To do this we introduce a force $\mathbf{F}_x$ acting on each particle $i$ in the nano-porous material.In the case of a gravitational field, we can formulate Darcy's law for flow as\begin{align}    \mathbf{U} = \frac{k}{\mu}\left( \nabla P - \rho \mathbf{a} \right).\end{align}Here $\mathbf{U}$ is the _Darcy velocity_, $k$ the permeabillity, $\mu$ the fluid viscosity, $P$ the externally applied pressure, $\rho$ the density and $\mathbf{a}$ the acceleration from the applied force. In our case $\nabla P = 0$ as we have no externally applied pressure. Looking at Newton's second law for a single particle we have\begin{align}    \mathbf{F}_x = m_i \mathbf{a}_i.\end{align}For all particles of the same mass with the same applied force this yields the total applied for to the system\begin{align}    N\mathbf{F}_x = N m \mathbf{a},\end{align}where $N$ is the number of particles. Dividng by the volume of the system we find\begin{gather}    \frac{N}{V}\mathbf{F}_x = \frac{N m}{V} \mathbf{a}    \implies    n \mathbf{F}_x = \rho \mathbf{a},\end{gather}where $n$ is the number density and $\rho$ the density of the system. Thus we see that we can replace $\rho \mathbf{a}$ in Darcy's law with $n\mathbf{F}_x$ and we write\begin{align}    \mathbf{U} = \frac{k}{\mu}\left( \nabla P - n \mathbf{F}_x \right).\end{align}<jupyter_code>import os
import sys

sys.path.append("..")

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

from read_lammps_dump import read_dump
from read_lammps_log import read_log, get_temp_lognames

sns.set(color_codes=True)
%%writefile scripts/flow_profile.in

# 3d Lennard-Jones gas
units lj
dimension 3
# Periodic boundiaries
boundary p p p
atom_style atomic

variable seed equal 87287
variable sigma equal 3.405
variable b equal 5.72
variable reduced_density equal 4/((${b}/${sigma})^3)
variable temperature equal 1.5
variable num_box equal 20
variable radius_min equal 20/${sigma}
variable radius_max equal 30/${sigma}

# Set fcc lattice with specified density
lattice fcc ${reduced_density}
region simbox block 0 ${num_box} 0 ${num_box} 0 ${num_box}
create_box 2 simbox
create_atoms 1 box

variable box_len equal lx

mass * 1.0
velocity all create ${temperature} ${seed}

pair_style lj/cut 3.0
pair_coeff * * 1.0 1.0

fix 1 all nvt temp ${temperature} ${temperature} 0.5

thermo 100
run 1000

# Use a radius of 2 nm (equal to 20 Å)
variable radius equal 20/${sigma}
variable centre_x equal lx/2
variable centre_y equal ly/2
variable centre_z equal lz/2
region cylinder_reg cylinder x ${centre_y} ${centre_z} ${radius} EDGE EDGE units box

# Delete half of the atoms in the region
delete_atoms porosity cylinder_reg 0.5 ${seed}

group pore_group region cylinder_reg
set group pore_group type 2

# Create another group from all the remaining atoms except for the ones in the cylinder
group frozen subtract all pore_group
# Set the velocity of the atoms outside the cylinder to zero
velocity frozen set 0 0 0
# Avoid integrating all the particles
unfix 1

# Run pore_group using NVE
fix 1 pore_group nve

dump 1 all custom 10 dat/flow_profile.lammpstrj id type x y z vx vy vz

# Compute the porosity
variable porosity equal count(pore_group)/count(all)
print ${porosity} file dat/flow_profile_porosity.dat

reset_timestep 0
compute pore_temp pore_group temp
compute msd pore_group msd
variable time equal dt*step

thermo_style custom step v_time c_pore_temp etotal press c_msd[4]

fix force pore_group addforce 0.1 0.0 0.0

thermo 10
log dat/flow_profile.log

run 2000
!export OMP_NUM_THREADS=4 && mpirun -np 4 lmp -in scripts/flow_profile.in
df, num_atoms, bounds = read_dump("dat/flow_profile.lammpstrj")
centre_bounds = [np.sum(bound) / 2 for bound in bounds]
df.head()
df_flowing = df[df.type == 2]
df_flowing.head()
df_flowing.x -= centre_bounds[0]
df_flowing.y -= centre_bounds[1]
df_flowing.z -= centre_bounds[2]
df_flowing.loc["r"] = np.sqrt(df_flowing["y"] ** 2 + df_flowing["z"] ** 2)
df_flowing.head()
r = df_flowing[df_flowing.timestep == np.max(df_flowing.timestep)].r

hist, bins = np.histogram(
    np.abs(df_flowing[df_flowing.timestep == np.max(df_flowing.timestep)].vx),
    np.linspace(0, np.max(r), 101)
)
fig = plt.figure(figsize=(14, 10))

plt.hist(hist, bins)
plt.show()
fig = plt.figure(figsize=(14, 10))

plt.scatter(
    df_flowing[df_flowing.timestep == np.max(df_flowing.timestep)].r,
    np.abs(df_flowing[df_flowing.timestep == np.max(df_flowing.timestep)].vx)
)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10492227979274611, 'fraction_numerical': 0.022020725388601035, 'mean_word_length': 4.712700369913687, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '63463', 'license': 'MIT', 'n_tokens_mistral': 1718, 'n_tokens_neox': 1586, 'n_words': 557, 'path': 'project-2/flow-in-a-nano-porous-material.ipynb', 'repo_name': 'Schoyen/FYS4460', 'size': '4632'}