Una grieta atravesando un cristal

Publicado el 17 de octubre de 2026

Una grieta no puede correr más deprisa que las ondas que avisan de que está llegando. Ese límite —la velocidad de Rayleigh— sale de las constantes elásticas del material, y en un cristal bidimensional de Lennard-Jones las constantes elásticas se pueden escribir en una línea. Así que aquí no hay nada que ajustar: la referencia es una fórmula cerrada y la simulación tiene que caer donde la teoría dice. Esto es ese cálculo completo.

Nota de notación. En unidades de Lennard-Jones hay dos colisiones de símbolos inevitables: $ϵ$ es la energía del potencial y $\epsilon$ la deformación; $\sigma$ es la distancia del potencial y ${\sigma }_{yy}$ la tensión. Las dos son convenio estándar y no merece la pena inventarse otro, así que aquí la deformación y la tensión van siempre con subíndice y los símbolos sueltos son las unidades.

Por qué se rompen los cristales, y por qué la teoría continua no basta

La mecánica de la fractura clásica responde a cuándo se rompe algo con un balance de energía. Una grieta avanza si al avanzar libera más energía elástica de la que cuesta crear las dos superficies nuevas. Eso es Griffith, y es una idea de 1921 que funciona asombrosamente bien: para una entalla de longitud $c$ en el borde de una lámina,

donde $\gamma$ es la energía de superficie, $E$ el módulo de Young y ese 1,12 el factor geométrico de una grieta que sale del borde en lugar de estar en el interior. Con los números de este cristal —que se validan más abajo— la deformación crítica sale ${\epsilon }_c\approx 0,0088$.

El problema de esa fórmula es que no contiene ni un átomo. Trata el sólido como un medio continuo, la grieta como un corte matemático de espesor cero y la superficie como una energía por unidad de longitud. Hay tres cosas que se le escapan, y las tres aparecen en cuanto se miran los átomos:

1. La punta no es un punto. La solución elástica da una tensión que diverge como $1/\sqrt{r}$ al acercarse a la punta. Eso no puede ser: a un par de distancias interatómicas de la punta la tensión es finita, porque los enlaces se estiran hasta que se rompen y no más. La zona donde la teoría continua deja de valer mide unos pocos nanómetros —en esta simulación, unos pocos σ—, y es justo la zona que decide si la grieta avanza.
2. La red atrapa la grieta. Avanzar una fila de átomos cuesta pasar una barrera, igual que mover una dislocación cuesta el Peierls. Griffith dice que por encima de ${\epsilon }_c$ la grieta debe correr; la red dice que hasta cierta sobrecarga la punta se queda enganchada entre dos filas. En esta lámina eso ocurre de verdad: por debajo de ${\epsilon }_{yy}\approx 0,02$ —más del doble de la crítica de Griffith— la grieta no arranca, y cuando arranca se para varias veces por el camino. Se ve a simple vista en la gráfica de posición frente a tiempo.
3. Hay un techo de velocidad, y no es el del sonido. Una grieta que avanza tiene que mover el campo de tensión que la rodea, y ese campo se reorganiza a la velocidad de las ondas elásticas. La condición de que el balance de energía siga cerrando impone que la punta no supere la velocidad de las ondas de Rayleigh $c_R$, que son las que viajan por una superficie libre: exactamente la superficie que la grieta va creando. Es un límite teórico limpio, y medirlo es la forma más directa de comprobar que la simulación está haciendo mecánica y no decoración.

Este montaje existe para mirar esas tres cosas con átomos contados, y a la vez para tener una referencia cerrada contra la que medirse. Es la mejor combinación posible para un cálculo de validación: un sistema lo bastante simple como para resolverlo a mano y lo bastante real como para romperse.

El montaje

Una lámina rectangular de red triangular, con una entalla abierta desde el borde izquierdo quitando una fila de átomos, y tracción perpendicular a la grieta (modo I):

Cuatro decisiones del montaje que no son obvias y que condicionan el resultado:

El radio de corte no es un detalle de eficiencia, es el modelo. En la red triangular relajada los seis primeros vecinos están a 1,122 σ y los segundos a 1,943 σ. Cortar en 1,5 σ deja dentro los primeros y fuera los segundos, y eso convierte el cristal en un problema resoluble: con solo primeros vecinos, el equilibrio exige ${\phi }^{\prime }(a)=0$ y todo lo demás sale en cerrado. Con un corte de 2,5 σ —el habitual en Lennard-Jones— el cristal sería más realista y la validación dejaría de ser exacta.

Se carga de golpe y no se sigue cargando. Es el protocolo clásico de fractura en MD: se aplica una deformación afín por encima de la crítica, se congelan las filas de los extremos para mantenerla y se deja correr. La grieta avanza con la energía elástica ya almacenada, que es lo que le da una velocidad estacionaria limpia en lugar de una que crece mientras se tira. El precio es que la energía se agota: por eso la grieta se frena antes de llegar al final de la lámina.

La carga es 4,5 veces la crítica de Griffith. No es generosidad: por debajo de unas dos veces y media la punta se queda atrapada en la red y no pasa nada en 240 τ. Cargar a ${\epsilon }_{yy}=0,040$ garantiza que la grieta corra, a cambio de que el resultado no diga nada fino sobre el umbral.

En $x$ no hay periodicidad. Con el eje largo periódico, la entalla se cerraría sobre sí misma por el borde y tendría dos puntas: saldrían dos grietas corriendo una hacia la otra y chocarían en medio. Los bordes libres en $x$ cuestan que haya superficie donde no debería haberla, pero dejan una sola punta, que es lo que se quiere medir.

Antes de nada, la validación

En el artículo del nanohilo de cobre la validación consistía en comprobar que un potencial ajustado reproducía las tablas a las que se ajustó, que es un filtro necesario pero flojo. Aquí es distinto y mejor: las referencias son fórmulas, no datos experimentales.

Con solo primeros vecinos, el parámetro de red sale de ${\phi }^{\prime }(a)=0$, es decir $a=2^{1/6}\sigma$, y la fórmula de Born para potenciales de pares en la red triangular se escribe entera:

lammps/in.elastico mide esas mismas magnitudes en LAMMPS, relajando la caja y deformándola por diferencias centradas, y python/valida.py las compara:

Esto no es «error del 1 %»: es coincidencia a nivel de redondeo. Y como las referencias son fórmulas y no medidas, lo que la tabla certifica es exactamente lo que interesa antes de simular: que la cadena de cálculo está bien montada. El fichero correcto, la red bien construida, las deformaciones del tamaño adecuado, las unidades sin un factor perdido.

Dos cosas que merecen un comentario:

La energía de superficie sale exacta, y tiene su gracia. Con solo primeros vecinos y ${\phi }^{\prime }(a)=0$, cada enlace está ya en su mínimo, así que al cortar el cristal los átomos del borde no sienten fuerza neta: la superficie no se relaja. El valor analítico —contar los enlaces que se rompen— es el valor correcto, y lo que queda de diferencia (0,000015 %) es la relajación residual que LAMMPS sí permite. En un cristal con más vecinos esto no pasaría, habría que medirla, y Griffith heredaría esa incertidumbre.

De la tabla sale todo lo demás sin tocar nada. El coeficiente de Poisson $\nu =C_{12}/C_{11}=1/3$ exacto; la isotropía elástica, que se comprueba con $C_{66}=(C_{11}-C_{12})/2$ y cuadra a $8\cdot 10^{-9}$; el módulo de Young bidimensional $E=(C_{11}^2-C_{12}^2)/C_{11}=65,987ϵ/{\sigma }^2$; y de ahí las velocidades de onda, que son las que de verdad hacen falta:

todas en σ/τ, con $\rho =0,916486$. Esa última, $c_R=4,840\sigma /\tau$, es la referencia contra la que se mide la grieta. No se ha ajustado nada para obtenerla: sale de una segunda derivada del potencial y de un poco de álgebra.

variable  a0     equal 2.0^(1.0/6.0)
variable  rho    equal 2.0/(sqrt(3.0)*2.0^(1.0/3.0))
lattice         custom ${rho} a1 1.0 0.0 0.0 a2 0.0 1.7320508075688772 0.0 &
                a3 0.0 0.0 1.0 basis 0.0 0.0 0.0 basis 0.5 0.5 0.0
pair_style      lj/cut 1.5
pair_coeff      1 1 1.0 1.0
pair_modify     shift yes

fix             relaja all box/relax x 0.0 y 0.0
min_style       cg
minimize        0.0 1.0e-14 100000 1000000
unfix           relaja
minimize        0.0 1.0e-14 100000 1000000

change_box      all x final 0.0 $((1.0+v_eps)*v_lx0) remap units box
minimize        0.0 1.0e-14 100000 1000000
variable  pxxp  equal pxx
variable  pxxP  equal ${pxxp}
change_box      all x final 0.0 ${lx0} remap units box
...
variable  C11   equal -(${pxxP}-${pxxM})/(2.0*${eps})

variable  epot0 equal pe
variable  e0    equal ${epot0}
change_box      all boundary p s p
minimize        0.0 1.0e-14 100000 1000000
variable  gamma equal (${epot1}-${e0})/(2.0*${lx0})

Qué sale

La animación es la simulación, no un vídeo: el navegador reconstruye los 63 fotogramas y dibuja los átomos en vivo a partir de un binario de 960 KiB. La barra de abajo recorre la trayectoria a mano.

tensión σ_yy (ε/σ²)03,2

Cómo leer los colores. Cada átomo va pintado según ${\sigma }_{yy}$, que es la componente que abre la grieta en modo I. Hay cuatro cosas que mirar, y las cuatro son física, no decoración:

• La lámina está entera en ámbar. Es la deformación aplicada de golpe: todo el cristal está cargado a la vez, a ${\epsilon }_{yy}=0,040$. (El primer fotograma, en t = 0, sale negro: la tensión que se pinta es una media de 20 pasos y en el instante cero todavía no hay ninguno. Es el único fotograma que no significa nada.)
• La entalla es una franja azul marino. Donde hay superficie libre no hay tensión que transmitir. Esa franja es la grieta inicial, y su longitud es el $c$ de la fórmula de Griffith.
• Delante de la punta hay un foco ámbar claro. Es la concentración de tensión, el $1/\sqrt{r}$ de la solución elástica, resuelto átomo a átomo. Es lo que rompe el enlace siguiente, y viaja con la punta.
• Detrás de la punta queda una estela oscura que se va ensanchando. La grieta descarga el material que va dejando atrás: ahí ya no hay nada que tire, porque los dos labios están separados. Esa estela es la energía elástica que se está consumiendo, y cuando el triángulo oscuro llega arriba y abajo de la lámina, se acabó el combustible.

La velocidad frente al límite de Rayleigh

El resultado cuantitativo es la trayectoria de la punta. Se localiza en cada volcado mirando la franja central de la lámina y avanzando desde el borde izquierdo mientras no haya átomos: donde vuelve a haberlos, está la punta.

Tres lecturas, por orden de lo que sostienen:

0,30 $c_R$ es un valor de libro para fractura frágil. La teoría dice que una grieta no puede superar $c_R$ —por encima, el balance de energía no se sostiene— y la experiencia dice que en la práctica se queda entre 0,2 y 0,5, con inestabilidades y ramificación cuando se acerca a la mitad. Que esta simulación caiga en 0,298 sin que nadie la haya empujado hacia ahí es la señal de que está haciendo mecánica de la fractura.

La escalera es el atrapamiento de red, en directo. La punta no avanza suave: se para en $x\approx 90\sigma$ durante 25 τ, corre, se vuelve a parar cerca de 225 σ. Cada descanso es la punta enganchada entre dos filas de átomos, esperando a que una fluctuación térmica o la onda que vuelve de un borde le dé el empujón que le falta. Esto es exactamente lo que la teoría continua no puede tener: en Griffith la grieta o avanza o no avanza, no duda. Es también la razón de que la velocidad se mida sobre el tramo central y no sobre el recorrido completo: la media global mezclaría propagación con descansos y daría 1,07 σ/τ, un número que no corresponde a ningún régimen.

El frenazo del final no es un resultado, es el protocolo. A partir de t ≈ 210 τ la punta se queda en 250 σ y ya no se mueve. No es que el cristal aguante más ahí: es que la energía elástica almacenada se ha consumido, que es lo que tiene que pasar cuando se carga de golpe y no se sigue cargando. Dónde se para depende de cuánto se cargó, no de una propiedad del material.

boundary        s s p
region          caja block 0 ${nx} 0 ${ny} -0.5 0.5 units lattice
create_atoms    1 box

variable  ymed   equal 0.5*(ylo+yhi)
variable  xent   equal xlo+${entalla}*lx
region    corte  block INF ${xent} $(v_ymed-0.35) $(v_ymed+0.35) INF INF units box
delete_atoms region corte compress yes

variable  ly0v   equal ly
variable  ly0    equal ${ly0v}
change_box      all y final $(ylo) $(ylo+(1.0+v_eps)*v_ly0) remap units box

group     borde  union sup inf
group     movil  subtract all borde
velocity  borde set 0.0 0.0 0.0
fix       congela borde setforce 0.0 0.0 0.0
fix       integra movil nve
fix       termo movil langevin ${temp} ${temp} 1.0 ${semi} zero yes
fix       dos all enforce2d

compute   sat  all stress/atom NULL
variable  Aat  equal sqrt(3)/2*${a0}^2
fix       prom all ave/atom 1 20 20 c_sat[2]
dump      d1 all custom ${nvolc} ${sal}.lammpstrj id x y f_prom

El campo de tensión, y por qué hay que suavizarlo dos veces

La tensión de un solo átomo fluctúa muchísimo más que el campo macroscópico: es una magnitud instantánea de un objeto que está vibrando. Pintada en crudo, la lámina es confeti y no se ve nada. Se promedia dos veces, y las dos veces por razones distintas:

• En el tiempo, con fix ave/atom 1 20 20 dentro de LAMMPS: cada valor volcado es la media de 20 pasos. Quita el ruido térmico de alta frecuencia, que no es información sobre el estado de tensión.
• En el espacio, con una gaussiana de 1,6 σ en python/campo.py, hecha con una rejilla para que sea $O(N)$. Esto no es cosmética: el campo de tensión del que habla la mecánica del continuo es un promedio sobre un volumen que contiene muchos átomos. Suavizar es construir la magnitud que se quiere comparar con la teoría, no disimular la que hay.

El mapa de color va de azul marino a ámbar pasando por arena, interpolando en espacio Lab con luminosidad creciente, de modo que la imagen se sigue leyendo en escala de grises o con daltonismo; los extremos son los dos colores de la marca. Y el rango se satura en ${\sigma }_{yy}=3,2$, no en el máximo: la singularidad de la punta es tan aguda que, escalando al máximo, todo lo demás queda plano. Es una decisión de visualización que cambia lo que se ve, así que va dicha: el ámbar más claro significa «3,2 ε/σ² o más».

Qué no demuestra esto

Cuatro límites, por orden de importancia:

La comparación con $c_R$ es con un límite, no con un valor. Esta es la más importante y la más fácil de vender de más. Las constantes elásticas y $\gamma$ sí son comparaciones exactas contra fórmulas cerradas: ahí se puede decir «coincide». La velocidad de la grieta no: la teoría dice $v<c_R$ y la simulación da $0,298c_R$. Eso es compatible con la teoría, no una predicción acertada. Una comprobación de verdad del régimen de velocidades pediría barrer la carga aplicada y ver que $v/c_R$ crece hacia el límite sin cruzarlo, que es lo que dice la teoría dinámica de la fractura y lo que este montaje no hace.

Griffith aquí solo da el orden de magnitud. Se carga a 4,5 veces la crítica porque por debajo de unas 2,5 veces la grieta no arranca. Es decir: el umbral real de esta lámina está bastante por encima del de Griffith, y la diferencia es el atrapamiento de red. La fórmula sigue siendo útil —dice en qué orden de magnitud buscar la carga— pero no es una predicción cuantitativa en un sistema de este tamaño, y presentarla como tal sería hacer trampa.

El cristal es bidimensional y el potencial es de juguete. Dos dimensiones no son tres: no hay dislocaciones de tornillo, no hay planos de deslizamiento que compitan, y la superficie nueva es una línea. Y un Lennard-Jones cortado entre vecinos no es ningún material: es el sistema modelo más simple que todavía se rompe. La ventaja es justamente esa —permite referencias exactas—, pero lo que se mide aquí son propiedades de este modelo. Lo que es transferible es el método y el mecanismo, no los números.

Una ejecución, una semilla. La grieta arranca y se para por sucesos individuales, y cuáles son depende de las velocidades iniciales. Para dar la velocidad con una barra de error —y para decir dónde está de verdad el umbral de arranque— habría que repetir con varias semillas y varias cargas. Son dos barridos que este montaje permite (EPS es una variable de entorno) y que este artículo no hace.

Nada de esto quita valor al cálculo; lo encuadra. Lo que sí sostiene este trabajo es la cadena entera: un cristal cuyas propiedades elásticas se conocen en cerrado y se reproducen a nivel de redondeo, una grieta que corre a una fracción sensata del límite teórico, y un campo de tensión que se puede mirar átomo a átomo y se comporta como la solución elástica dice.

Cómo reproducirlo

Todo está en github.com/djoroya/sim-visuales, en md/grieta-2d/. Un solo script:

./run.sh              # validación + simulación + salidas
./run.sh validacion   # solo las constantes elásticas y γ
./run.sh simulacion   # solo el ensayo de fractura (unos 3 minutos)
./run.sh salidas      # vídeo, imágenes y datos del visor

Hace falta LAMMPS sin paquetes especiales —lj/cut y enforce2d van en el núcleo—, ffmpeg, y Python 3 con numpy, matplotlib y Pillow; el script crea el entorno virtual la primera vez. Las variables útiles son NX y NY (celdas de la lámina), EPS (la deformación aplicada), TFIN y HILOS.

El volcado completo son 283 MB y se queda en datos/, fuera de git. Lo que se versiona es lo reproducible: las dos entradas de LAMMPS, los scripts de Python y las salidas pequeñas.

Etiquetas: dinámica molecular · LAMMPS · fractura · validación