Investigación / Fractura y daño

Fractura por campo de fase de una placa entallada

A partir de Nathan Shauer, “Less than 500 lines self-contained Python finite element implementation of the phase-field method for fracture mechanics”, Advances in Engineering Software (2025), doi:10.1016/j.advengsoft.2025.104013.

Frente al código del autor
30 de 30 controles aprobados· ver los controles
Frente a las figuras del artículo
No comparado: no tuvimos acceso al texto completo

Investigación, no una herramienta de diseño. Nada de esto es un resultado de ingeniería.

Una placa cuadrada con una entalla desde su borde izquierdo hasta el centro se tracciona desde arriba. El modelo de campo de fase AT2 representa la grieta como un campo de daño d entre 0 (intacto) y 1 (roto), distribuido sobre una longitud ℓ. El panel del experimento ejecuta la implementación Rust de Stabileo en su visor; la curva naranja es lo que calculó el código Python del autor.

Qué muestra el artículo

El artículo presenta un código autónomo de elementos finitos en Python, de menos de 500 líneas, para la fractura frágil con el método de campo de fase, y lo aplica a tres ejemplos: una barra a tracción comparada con la solución unidimensional AT2, una placa con entalla en un borde sometida a tracción (esta página) y un ensayo de flexión en tres puntos. Leímos todo el código (no tuvimos acceso al texto del artículo) y todo lo que sigue sobre el método proviene de él.

En la reproducción el daño crece suavemente por delante de la entalla mientras crece linealmente hasta ; en el siguiente incremento , después de , y la fuerza cae casi a cero. muestra la grieta abierta en la placa deformada.

El modelo

Desplazamiento u y el daño d son ambos bilineales en cuadriláteros, integrados con 2×2 puntos de Gauss; deformación plana, pequeñas deformaciones, carga cuasiestática.

  • Elasticidad con d fijo: div(g(d) C : ε(u)) = 0 con g(d) = (1 − d)², aplicado a toda la tensión.
  • Daño (AT2) con valor fijo de u: (Gc/ℓ + 2H) d − Gc ℓ Δd = 2H.
  • Energía impulsora: H = ψ+ = ½ Σ σi εi sobre las deformaciones principales positivas, con σi las tensiones principales sin daño: la elasticidad se degrada por completo mientras que el daño es impulsado solo por la tracción, una formulación «híbrida».
  • Solución alternada por paso de carga: resolver u con el d actual, luego d con el nuevo u, hasta que el residuo de elasticidad ‖K(d)u − f‖ sea menor que 10−5.
  • Irreversibilidad: ninguna en el código de referencia; el daño se recalcula a partir de la energía actual y no está acotado por 1. para conservarlo durante la descarga (no está en la referencia).

Las unidades son kN y mm. La tensión máxima homogénea AT2 σc = √(27 E Gc / 256 ℓ) ≈ 2.82 kN/mm² aquí, por lo que ℓ actúa tanto como una resistencia como una longitud de regularización: y la carga máxima baja. La malla también importa: la banda de la grieta necesita elementos varias veces más pequeños que ℓ; .

Cómo lo reprodujimos

  1. El código del autor (commit 61037c6, licencia MIT) se ejecutó sin modificaciones para sus tres ejemplos en Python 3.12.15 con numpy 2.1.3 y scipy 1.14.1; la fuerza, el número de iteraciones, el residuo y el rango de daño de cada paso, y los campos nodales en pasos seleccionados, se guardaron con un manifiesto de versiones y sumas de comprobación.
  2. La implementación Rust regenera las mismas mallas (idénticas bit a bit), aplica las mismas condiciones de contorno, bucle y criterio de parada, y resuelve con la factorización de Cholesky dispersa del motor de Stabileo.
  3. 30 controles con nombre comparan ambas implementaciones, cada uno con una tolerancia declarada, en integración continua. El paquete del navegador y la reproducción de esta página se verifican por separado frente a la compilación nativa.

Comparación con la referencia

30 controles de validación ejecutados en CI, cada uno con su tolerancia
ControlObservadoToleranciaResultado
ex1: mesh identical to the reference (node and element fingerprints)igual—✓ aprobado
ex1: same number of load stepsigual—✓ aprobado
ex1: same number of staggered iterations at every stepigual—✓ aprobado
ex1: maximum damage at every step8,3 × 10−121 × 10−9✓ aprobado
ex1: stress at the bar centre at every step1,1 × 10−121 × 10−9✓ aprobado
ex1: nodal damage at step 504,5 × 10−101 × 10−8✓ aprobado
ex1: nodal damage at step 733,8 × 10−101 × 10−8✓ aprobado
ex1: nodal displacement at step 501,6 × 10−101 × 10−8✓ aprobado
ex1: centre stress vs 1D AT2 homogeneous solution while max damage < 0.257,3 × 10−51 × 10−3✓ aprobado
ex2: mesh identical to the reference (node and element fingerprints)igual—✓ aprobado
ex2: same number of load stepsigual—✓ aprobado
ex2: same number of staggered iterations at every stepigual—✓ aprobado
ex2: peak force at the same load stepigual—✓ aprobado
ex2: peak force1,2 × 10−121 × 10−9✓ aprobado
ex2: force at every step up to the peak1,6 × 10−121 × 10−9✓ aprobado
ex2: maximum damage at every step up to the peak4,2 × 10−121 × 10−9✓ aprobado
ex2: force at every step, including fracture (comparison tolerance)2,2 × 10−71 × 10−5✓ aprobado
ex2: nodal damage at the peak step 565 × 10−111 × 10−8✓ aprobado
ex2: nodal displacement at the peak step 568,5 × 10−111 × 10−8✓ aprobado
ex2: nodal damage at step 57, reference or its mirror image2,9 × 10−51 × 10−4✓ aprobado
ex2: nodal damage at step 78, reference or its mirror image3,1 × 10−101 × 10−8✓ aprobado
ex3: mesh identical to the reference (node and element fingerprints)igual—✓ aprobado
ex3: same number of load stepsigual—✓ aprobado
ex3: same number of staggered iterations at every stepigual—✓ aprobado
ex3: peak force at the same load stepigual—✓ aprobado
ex3: force at every step up to the peak4 × 10−111 × 10−9✓ aprobado
ex3: force at every step, including fracture (comparison tolerance)1 × 10−61 × 10−4✓ aprobado
ex3: nodal damage at the peak step 566,9 × 10−111 × 10−8✓ aprobado
ex3: nodal displacement at the peak step 561,1 × 10−101 × 10−8✓ aprobado
ex3: nodal damage at the last step 98, reference or its mirror image4,1 × 10−101 × 10−8✓ aprobado

Un hallazgo: la placa es simétrica respecto al plano de la entalla, y la grieta avanza por la fila de elementos justo encima en un código y justo debajo en el otro. La fila elegida depende del redondeo en coma flotante; frente a la imagen especular los campos de daño finales difieren en 3.085970368132962 × 10−10. Los umbrales de la tabla son tolerancias de comparación, no cotas matemáticas del error. Los campos de desplazamiento se comparan en el ejemplo 1 paso 50, el ejemplo 2 paso 56 y el ejemplo 3 paso 56.

Tabla de pasos (79 pasos)
pasou, mmF, kNiteracionesresiduomax dF de referencia, kN
00,000100,0134528712,8 × 10−80,00010,01345287
10,000200,0269047511,7 × 10−70,00040,02690475
20,000300,0403544314,3 × 10−70,00090,04035443
30,000400,053800717,9 × 10−70,00160,0538007
40,000500,0672423611,3 × 10−60,00240,06724236
50,000600,0806781911,9 × 10−60,00350,08067819
60,000700,0941069812,6 × 10−60,00480,09410698
70,000800,107527513,4 × 10−60,00630,1075275
80,000900,120938614,4 × 10−60,00790,1209386
90,001000,13433915,5 × 10−60,00980,134339
100,001100,147727616,7 × 10−60,01190,1477276
110,001200,16110318 × 10−60,01420,161103
120,001300,174464219,4 × 10−60,01660,1744642
130,001400,187788723,5 × 10−70,01940,1877887
140,001500,201114324,5 × 10−70,02230,2011143
150,001600,214421825,9 × 10−70,02550,2144218
160,001700,227709827,6 × 10−70,02890,2277098
170,001800,240977129,6 × 10−70,03240,2409771
180,001900,254222621,2 × 10−60,03620,2542226
190,002000,26744521,5 × 10−60,04030,267445
200,002100,280642921,8 × 10−60,04460,2806429
210,002200,293815322,2 × 10−60,04910,2938153
220,002300,306960822,7 × 10−60,05390,3069608
230,002400,320078223,2 × 10−60,05890,3200782
240,002500,333166223,8 × 10−60,06410,3331662
250,002600,346223524,5 × 10−60,06970,3462235
260,002700,359248925,3 × 10−60,07550,3592489
270,002800,37224126,1 × 10−60,08160,372241
280,002900,385198627,1 × 10−60,08800,3851986
290,003000,398120328,3 × 10−60,09470,3981203
300,003100,411004829,6 × 10−60,10170,4110048
310,003200,423841132 × 10−60,10920,4238411
320,003300,436645632,4 × 10−60,11680,4366456
330,003400,449408632,9 × 10−60,12490,4494086
340,003500,462128633,5 × 10−60,13330,4621286
350,003600,474804134,3 × 10−60,14210,4748041
360,003700,487433735,1 × 10−60,15130,4874337
370,003800,500015636,2 × 10−60,16090,5000156
380,003900,512548337,4 × 10−60,17100,5125483
390,004000,5250338,9 × 10−60,18160,52503
400,004100,53745243,3 × 10−60,19290,537452
410,004200,549824944,1 × 10−60,20460,5498249
420,004300,56214145,1 × 10−60,21700,562141
430,004400,57439846,5 × 10−60,22990,574398
440,004500,586593548,2 × 10−60,24370,5865935
450,004600,598718754,2 × 10−60,25840,5987187
460,004700,610781455,5 × 10−60,27390,6107814
470,004800,622773357,4 × 10−60,29030,6227733
480,004900,634690459,9 × 10−60,30800,6346904
490,005000,646520566,8 × 10−60,32720,6465205
500,005100,6582769,7 × 10−60,34790,65827
510,005200,669918178 × 10−60,37080,6699181
520,005300,681459687,6 × 10−60,39620,6814596
530,005400,692876598,5 × 10−60,42500,6928765
540,005500,704134118,1 × 10−60,45880,704134
550,005600,7151725148,7 × 10−60,50060,7151725
560,005700,7257948229,3 × 10−60,56060,7257948
570,005800,00047543893088,3 × 10−61,00810,0004752194
580,005900,000454193924,6 × 10−71,00810,0004541974
590,006000,000433467822,5 × 10−71,00830,0004334679
600,006100,000413999222,3 × 10−71,00850,0004139992
610,006200,000395713222,1 × 10−71,00870,0003957132
620,006300,000378519822 × 10−71,00880,0003785198
630,006400,000362336521,9 × 10−71,00890,0003623365
640,006500,000367766319,6 × 10−61,00910,0003677663
650,006600,000352428319,2 × 10−61,00920,0003524283
660,006700,000337749318,9 × 10−61,00930,0003377493
670,006800,000323894718,5 × 10−61,00940,0003238947
680,006900,000310807618,2 × 10−61,00940,0003108076
690,007000,00029843417,9 × 10−61,00950,000298434
700,007100,000286724817,6 × 10−61,00950,0002867248
710,007200,000275635117,4 × 10−61,00950,0002756351
720,007300,000265123317,1 × 10−61,00960,0002651233
730,007400,000255151616,9 × 10−61,00960,0002551516
740,007500,000245684916,7 × 10−61,00960,0002456849
750,007600,00023669116,4 × 10−61,00960,000236691
760,007700,0002281416,2 × 10−61,00960,00022814
770,007800,000220004516 × 10−61,00960,0002200045
780,007900,00021225915,8 × 10−61,00950,000212259

Comparación con el artículo

No realizada: no tuvimos acceso al texto completo del artículo, por lo que no se comparó ningún número publicado en él. Las figuras aquí provienen de nuestras propias ejecuciones.

Limitaciones

  • La concordancia con un código de referencia no es una validación frente a experimentos, y nada de esto es una certificación de ingeniería.
  • La malla de referencia refina solo una caja de 0.04 × 0.04 mm alrededor de la punta de la entalla; fuera de ella la grieta cruza elementos más anchos que 3ℓ.
  • Sin restricciones de cota: el daño puede superar ligeramente 1 (la referencia alcanza 1.0096 aquí y 1.0021897551918642 en el ejemplo 1) y caer por debajo de 0. Los resultados de referencia Python registrados tienen daños nodales mínimos de 3.076336950017797e-08 (ejemplo 1), -0.0031224948220514573 (ejemplo 2) y -0.06024191772700118 (ejemplo 3). Esto también ocurre en la discretización AT2 sin restricciones del código de referencia. El daño negativo no es físico: (1-d)² en el mínimo del ejemplo 3 es aproximadamente 1.124, por encima de la rigidez intacta. El ensamblaje usa daño interpolado en puntos de cuadratura, así que esto no es un error de 12.4% en la rigidez ensamblada. El solver no recorta el daño; el campo nodal cuantizado de la reproducción y la escala de colores lo recortan para su visualización.
  • El esquema alternado no está acelerado; durante el crecimiento inestable de la grieta puede requerir cientos de iteraciones en un paso.
  • Las ejecuciones en el navegador están limitadas en tamaño de malla y rango de parámetros.

Código y referencias

  • Implementación, pruebas y resultados de referencia: shauer-2025 (generadores de malla y bucle adaptados de la referencia MIT, con atribución en el código fuente).
  • Código de referencia: github.com/nathanshauer/phasefield-jr-py, MIT, © 2025 Nathan Shauer.

Citar esto

Stabileo Research, “Fractura por campo de fase de una placa entallada”, reproducción de Nathan Shauer (2025), doi:10.1016/j.advengsoft.2025.104013. Implementación: research/fracture/shauer-2025 en el commit 0ce459e9eb75. Motor: el motor de Stabileo en el mismo commit. Se prevén enlaces permanentes a esta compilación exacta.