La principal diferencia entre el enfoque de implementación directa vs. el enfoque de Multifísica discreta (DMP) radica en la complejidad del primer enfoque en contraste con la simplicidad de adaptación del enfoque DMP. En el primer caso, el de implementación de estilo de par directamente en LAMMPS, hay que generar toda una nueva clase que deberá ajustarse a algún modelo de comportamiento de fluido viscoelástico, por ejemplo mediante la selección de uno de los modelos estándar presentados en 2, y la genereación de una nueva clase de estilo de par. Este proceso conlleva redefinir todos los métodos requeridos por el nuevo estilo de par y, por tanto, resultar en un trabajo que puede ser bastante intensivo o dispendioso. Por otro lado, mediante el enfoque DMP basta con hacer una selección cuidadosa de los estilos de par a ser utilizados de manera concurrente, y generar una adecuada definición de caso en el script de caso de LAMMPS. En cualquier caso, independiente del enfoque de simulación empleado, es necesario llevar a cabo un proceso de selección de los parámetros de simulación de los pares de estilo seleccionados para llevara a cabo la simulación de flujo viscoelástico. A continuación se presentan algunos detalles específicos del tratamiento que debe hacerse en cada uno de estos enfoques.
Como caso de contexto de implementación de un modelo clásico viscoelástico se presentaran algunos detalles de la implementación directa del modelo de Kelvin. Para este fin se usaran los conceptos mostrados en la sección anterior (3) para generar un nuevo potencial de enlace disipativo que puede utilizarse para modelar sustancias viscoelásticos. Como se había indicado anteriormente, el modelo de Kelvin se usa para la modelación de materiales viscoelásticos mediante la suposición de un comportamiento equivalente al de un arreglo en paralelo de un componente tipo amortiguador puramente viscoso y de un componente tipo resorte puramente elástico conectados en paralelo, como se muestra en la .
Como se indicaba anteriormente, este arreglo de elementos disipativos y elasticos puede expresarse como: \[\label{eq:modeloKelvinV2} \sigma(t) = k \epsilon(t) + \eta \dot{\epsilon}\] donde el esfuerzo queda expresado como función de la deformación \(\epsilon\), la rata de deformanción \(\dot{\epsilon}\), y dos parámetros \(k\) y \(\eta\), representando la constante elástica y la constante disipativa viscosa, respectivamente.
Para llevar a cabo la implementación de un estilo de par, o de cualquier otra clase de hecho, es conveniente tomar como punto de partida una clase similar o asociada ya existente. Esta práctica es común en muchas aplicaciones de programación cientifica, con otro ejemplo extendido de esta práctica la programación de funciones o solucionadores personalizados usados en el programa , por mencionar un caso adicional. Para la implementación de este estilo de par personalizado se ha seleccionado como referencia la clase del estilo de par denominado , cuya estructura de implementación se presenta en , mostrado a continuación:
lst:bondHarmonicCppbondHarmonic.cpp
Como es costumbre, la nueva clase debe declararse e inicializarse mediante los archivos de encabezado, , y de código los que, en general, deberán guardarse en el directorio /src/MOLECULE de la carpeta de instalación de LAMMPS, y siguiendo una jerarquia como la mostrada en .
Todas las funciones serán las mismas que en el estilo de referencia. Sin embargo, en el nuevo estilo , tenemos que sustituir el texto “BondHarmonic" por un nuevo texto “BondKelvin", como puede verse en .
lst:bondKelvinCppbondKelvin.cpp
En comparación con el enlace de tipo armónico ahora se requiere de un nuevo parámetro, \(\eta\) o , desde el archivo de entrada. Por esta razón necesitamos modificar los métodos , , , , y .
Siguiendo el orden de inicialización de funciones (ver ), a continuación se muestran los códigos abreviados de estos métodos, en el orden presentado, para el estilo original y los códigos abreviados para el nuevo estilo .
lst:destructorBondHarmonicCppdestructorBondHarmonic.cpp
lst:destructorBondKelvinCppdestructorBondKelvin.cpp
El siguiente método a generar, siguiendo el orden presentado en , es . En el nuevo estilo de par se requiere la rata de deformación, \(\dot{\epsilon}\), que puede también considerarse como la velocidad de deformación. Para poder utilizar dicha variable dentro del nuevo estilo de pares necesitamos declarar e inicializar las velocidades de cada partícula, lo que constituye la modificación a ser hecha, en comparación con el estilo de par ármonico.
lst:computeBondHarmonicCppcomputeBondHarmonic.cpp
lst:computeBondKelvinCppcomputeBondKelvin.cpp
Dentro del bucle del método original, es necesario añadir un nuevo conjunto de líneas de código para calcular la fuerza del componente viscoso, así como las variaciones de energía. Este cáculo, que no se presenta en detalle, requiere de la estimación de velocidades de partícula y direcciones para cálculo de fuerzas del componente viscoso. Estas definiciones de variables se presentan en .
lst:computeVelsAndDirsBondKelvinCppcomputeVelsAndDirsBondKelvin.cpp
de manera que ahora es posible escribir las expresiones modificadas para el cálculo de la fuerza aplicada a cada par de partículas, como se presenta a continuación. Como en los caoss anteriores se presenta el código resumido para el estilo de par de referencia () y para el nuevo estilo de par ().
lst:computeForcesBondHarmonicCppcomputeForcesBondHarmonic.cpp
lst:computeForcesBondKelvinCppcomputeForcesBondKelvin.cpp
Siguiendo la línea de trabajo anterior, y dado que fue necesario introducir un nuevo parámetro en el estilo de par, es necesario modificar el método para considerar la necesidad de reserva de nueva memoria dinámica, lo que se presenta a continuación.
lst:allocateBondHarmonicCppallocateBondHarmonic.cpp
lst:allocateBondKelvinCppallocateBondKelvin.cpp
El coeficiente del componente de disipación viscosa \(\eta\) (o el análogo a un elemento amortiguador), debe ser suministrado por el usuario mediante línea de código en el archivo de configuración de caso de LAMMPS, por lo que es necesario también hacer ajustes enla definición del método , que se presenta en .
lst:coeffBondHarmonicCppcoeffBondHarmonic.cpp
lst:coeffBondKelvinCppcoeffBondKelvin.cpp
El estilo de par de referencia tiene también las funciones y que tienen que ser modificadas. Estas funciones de tipo escriben y leen un archivo de geometría que puede utilizarse como archivo de soporte en el archivo de entrada o de configuración de caso LAMMPS. El código de referencia y su respectiva modificación para estos métodos se presentan en
En el archivo encabezado del nuevo estilo de par se debe sustituir el texto “BondHarmonic" por un nuevo texto “BondKelvin", así como declarar un nuevo miembro protegido en la clase, el puntero a .
lst:bondHarmonicHbondHarmonic.h
lst:bondKelvinHbondKelvin.h
Finalmente, dado que el nuevo estilo debe invocarse desde el script de configuración de caso de LAMMPS, y ya que el nuevo estilo de par está completamente programado, es necesario compilarlo y luego invocarlo escribiendo las siguientes líneas de comando en el script de configuración de caso:
Mientras que la simulación de flujos viscoelásticos con implementación directa requirió la creación y ajuste de toda una nueva clase de estilo de par, con todo lo que implicó, el enfoque DMP se basa en llevar a cabo la simulación del flujo mediante acople de potenciales simples ya existentes, pero que emulen los componentes elástico y viscoso requeridos. En general, como se expondrá a continuación, se desea llevara a cabo una simple superposición aditiva de potenciales, de manera que el efecto final resultante será precisamente el de una sustancia con comportamiento elástico y viscoso simultaneamente.
Como se ha indicado en forma general, la viscoelasticidad puede implementarse directamente en SPH (por ejemplo, ). Sin embargo, de manera equivalente a la implementación directa de modelos convencionales de viscoelasticidad, esto implica reescribir la ecuación de movimiento para tener en cuenta un modelo específico de viscoelasticidad. En el caso de un código de partículas como LAMMPS, por ejemplo, esto implicaría reescribir grandes secciones del código dedicadas a SPH cada vez que introduzcamos un nuevo modelo de viscoelasticidad. El enfoque propuesto en este estudio modela la viscoelasticidad mediante la combinación de diferentes potenciales de partículas, que es un procedimiento estándar en los códigos de partículas, en lugar de reescribir las ecuaciones de movimiento. Para la parte viscosa, se adopta el enfoque SPH estándar descrito en . De forma análoga a la idea en la que se basan los modelos viscoelásticos básicos más utilizados anteriormente, se propone una técnica de modelación alternativa dentro de un marco basado en partículas, en el que los potenciales entre partículas que abordan las interacciones viscosas y elásticas por separado se mezclan de forma aditiva para imitar la respuesta dual disipativa y restauradora de las sustancias viscoelásticas. Con este fin, el soporte viscoso SPH se mezcla con un potencial que se asemeja a un comportamiento elástico restaurador, con el objetivo de obtener la respuesta dual característica de las sustancias viscoelásticas. En nuestro planteamiento, las propiedades viscosas vienen dadas al fluido por las ecuaciones SPH . Además, en principio, se busca utilizar potenciales equivalente a los resortes armónicos que se emplean en el Lattice Spring Model (LSM) para añadir elasticidad al material (véase ),
\[\label{eq:modelLSM} U_{\textrm{LSM}}(r) = \dfrac{1}{2} k (r-r_0)^2\]
que proporcionaría una fuerza \(F(r)=-\nabla U(r) = -k(r-r_0)\) entre dos partículas conectadas por el resorte con una rigidez \(k\), separadas una distancia \(r\), y con una distancia de equilibrio \(r_0\). La combinación demodelo SPH y el modelo tipo LSM conferirá propiedades viscoelásticas al material. En se empleó un enfoque similar para la modelación de sólidos viscoelásticos. Sin embargo, los fluidos no pueden modelarse del mismo modo porque los resortes limitan las partículas, que no pueden fluir como debería ocurrir en los fluidos. La solución que se propone en el presente trabajo es emplear una mezcla de dos potenciales para imitar parcialmente el LSM, pero con la ventaja de tener una distancia de desactivación o de corte a partir de la cual el potencial compuesto deja de actuar, permitiendo que las partículas fluyan libremente. Se eligió un potencial exclusivamente atractivo, propuesto por primera vez por , para proporcionar la parte atractiva del análogo del LSM, y definido como,
\[\label{eq:cosineSquaredDefinition} U_{\textrm{CS}}(r) = \left\{ \begin{array}{lll} -\epsCosSq & \quad r < \sigma \\ -\epsCosSq \cos{\left(\dfrac{\pi(r-\sigma)}{2(\rcCosSq-\sigma)}\right)}^{2} & \quad \sigma \leq r < \rcCosSq \\ 0 & \quad r \geq \rcCosSq \end{array} \right.\]
que, como se ilustra esquemáticamente en la , es un valor constante de por debajo de una distancia interparticular \(\sigma\), que aumenta proporcionalmente a \(r\) hasta que desaparece por encima de una distancia de corte , según la Ecuación . Esta enfoque permite tener en cuenta el hecho de que las partículas pueden dejar de sentir una fuerza de atracción cuando están muy separadas; así, las partículas de fluido no están constreñidas en una estructura reticular como en la ecuación , sino que son libres de fluir dentro del dominio.
Para evitar la superposición de partículas, en nuestro modelo la parte repulsiva se proporciona mediante otro potencial dado como, \[U_{\textrm{Soft}}(r) = \epsSoft \left[ 1 + \cos{\left(\frac{\pi r}{ \rcSoft }\right)} \right] \qquad r < \rcSoft\] donde es la magnitud y es la distancia de corte para este componente repulsivo del potencial adoptado, limitando así, una vez más, la distancia a la que existe repulsión para partículas muy alejadas. En el presente modelo propuesto utilizamos un potencial total añadiendo estos componentes con el objetivo de imitar el comportamiento elástico del modelo LSM, pero con la ventaja de que cuando la distancia entre las dos partículas está por encima del valor de corte, la fuerza se desactiva y la partícula es libre de fluir. En la se presenta una ilustración de los potenciales LSM, repulsivo, atractivo y total, mientras que en la se muestra una representación de las fuerzas producidas por ellos, donde se representan esquemáticamente los componentes repulsivo y atractivo de las fuerzas elásticas equivalentes.
Dado que la lógica del modelo propuesto es poder obtener un acoplamiento entre los dos estilos de pares que produzca una repulsión/atracción estable dentro de un rango cercano de cada partícula, este acoplamiento requiere una elección sensata de valores dada la posible interacción compleja que puede producirse por el número de parámetros en juego. Por ejemplo, una consideración importante, normalmente crítica en otros métodos basados en mallas sólo desde la perspectiva de la estabilidad numérica, es la resolución espacial. En nuestro caso, la resolución espacial está vinculada de algún modo a la distribución de las partículas, por lo que la estabilidad y eficacia del modelo depende también de la malla empleada para la disposición inicial de dichas partículas, así como de la longitud característica en la malla, o escala de malla \(\Delta_{L}\). Concretamente, con el fin de reducir el número de parámetros libres, el potencial repulsivo se definió en función de las características del potencial atractivo, de modo que un único conjunto de valores pudiera definir por completo el modelo acoplado elástico. Es importante mencionar que, dado que nuestro modelo depende fuertemente del espaciado entre partículas, los datos presentados a lo largo del presente trabajo como valores de referencia deben tomarse como una guía para establecer simulaciones de trabajo, más que como un rango restrictivo de condiciones de operación. En cualquier caso, algunos rangos y definiciones empleados en la caracterización de la componente elástica del presente modelo, determinados a través de varios experimentos numéricos y que se ha comprobado que proporcionan estabilidad numérica y consistencia con el comportamiento esperado, se presentan en la a modo de referencia. En esta tabla los acrónimos BCC y FCC significan Body-centred cubic (cúbica centrada en el cuerpo) y Face-centred cubic (cúbica centrada en la cara), respectivamente, que son algunos de los tipos más comunes de celdas regulares para simulaciones de dinámica molecular.
| BCC Lattice | FCC Lattice | |
|---|---|---|
| Distancia de activación - potencial atractivo | \(\sigma = \sqrt{3}\,\Delta_{L}/2\) | \(\sigma = \sqrt{2}\,\Delta_{L}/2\) |
| Distancia de corte - potencial atractivo | \(0.95\,\Delta_{L}\leq\rcCosSq\leq2.1\,\Delta_{L}\) | \(0.9\,\Delta_{L} \leq\rcCosSq\leq1.1\,\Delta_{L}\) |
| Distancia de corte - potencial repulsivo | \(\rcSoft \approx 1.05\,\rcCosSq\) | \(\rcSoft\approx 0.955\,\rcCosSq\) |
| Prefactor - potencial repulsivo | \(0.8\,\epsCosSq\leq\epsSoft\leq3.0\,\epsCosSq\) | \(1.25\,\epsCosSq\leq\epsSoft\leq4.0\,\epsCosSq\) |
Una vez definidos el tipo de estructura retícular a usar en la malla inicial (BCC o FCC) y los valores de los paraámetros de los potenciales atractivo y repulsivo del componente elástico, es necesario configurar un script de caso de simulación LAMMPS. Existen muchas posibilidades de estructura de tipo de archivo (). Un ejemplo de archivo de configuración de caso LAMMPS para ser usado en simulación de flujos viscoelásticos bajo el enfoque propuesto en este trabajo se presenta en el siguiente código().
La definición de los parámetros a usar para la simulación de casos de flujo viscoelástico mediante el enfoque DMP propuesto en el presente trabajo se puede encontrar entre las líneas y del . El resto del contenido del archivo de configuración es usado para definir geometría (líneas a ), parámetros de ejecución (líneas a ), dominio y coondiciones de frontera (líneas a ), propiedades de partícula (líneas a ) e instrucciones de posprocesamiento y ejecución (líneas a )