Una rotación rígida de 30° marcó -13% de deformación — pares conjugados tensión-deformación en láminas

En grandes deformaciones no se puede multiplicar cualquier tensión por cualquier deformación. Solo $S:\dot{E}$, $P:\dot{F}$ y $J\sigma:d$ dan el mismo valor.

Solo hubo rotación y el medidor de deformación marcó -13%#

Un elemento de lámina se giró 30° dentro de su propio plano. Sin estirarlo ni retorcerlo. Solo rotación rígida.

Aun así, la deformación ingenieril $\varepsilon_{xx} = \partial u_x / \partial x$ da $-0.134$ . Es un 13.4% de compresión. El elemento no se deformó en ningún punto y el medidor lee compresión.

Ese valor no es un descuido ni un error de discretización. $\cos 30° - 1 = -0.134$ , exactamente ese valor. La definición misma de la deformación infinitesimal confunde rotación con deformación.

Este texto trata de dónde se corta esa confusión y qué tensión hay que usar una vez cortada. Hay tres cosas por llevarse. Una medida de deformación inmune a la rotación, valores medidos de cuánto se desvía la energía cuando el par tensión-deformación está mal elegido, y la identidad de la ecuación que cierra la incógnita en el espesor de una lámina.

Con rigid rotation only activo, al barrer el ángulo solo se mueve la $\varepsilon$ roja. La $E$ verde no sale de cero hasta que sube lam.

theta30°lam1.00gam0.00

Press rigid rotation only and let it spin: the patch never changes shape, yet the red eps bars swing all the way across while the green E bars stay pinned at zero. At theta = 30° the small-strain gauge reads 0.0000 against 0.0000 — a gap of 0.0000 invented by the rotation alone. Now add lam and gam: E moves, and it keeps the same value at every theta.

La forma del elemento no cambia. La barra roja llega cerca de -0.134 a 30°, y la verde se queda en cero aunque el ángulo recorra todo el dial.

Lo que filtra la rotación no es $F$ sino $F^{T}F$ #

La deformación arranca con el gradiente de deformación (deformation gradient, la aplicación local de la configuración de referencia a la actual).

F_{iJ} = \frac{\partial x_i}{\partial X_J}

$x$ es la coordenada deformada y $X$ la no deformada. La descomposición polar (polar decomposition) la separa como $F = RU$ . $R$ es rotación y $U$ es estiramiento puro. El problema es que $R$ sigue vivo dentro del propio $F$ .

$\varepsilon = \tfrac{1}{2}(F + F^{T}) - I$ usa $F$ tal cual. Por eso $R$ se cuela. De ahí que bajo rotación pura resulte $\varepsilon_{xx} = \cos\theta - 1$ .

Pero al elevar $F$ al cuadrado la rotación desaparece.

C = F^{T}F = U^{T}R^{T}RU = U^{T}U

Como $R^{T}R = I$ , en $C$ solo queda $U$ . Ese $C$ es el tensor de Cauchy-Green por la derecha. Restarle la identidad y dividir por dos da la deformación de Green-Lagrange.

E = \tfrac{1}{2}\left(F^{T}F - I\right)

Bajo rotación rígida $F^{T}F = I$ , así que $E$ es exactamente cero. Esas dos líneas explican por qué la barra verde de la simulación no se mueve. La exigencia de que la física no cambie al cambiar el sistema de coordenadas se resolvió con un cambio de base en la transformación del tensor constitutivo; aquí se resuelve con la definición misma de la medida de deformación.

¿Por qué área está dividido el tensor de tensión?#

Si la deformación se llevó a la configuración de referencia, la tensión debe ir con ella. La tensión es "fuerza sobre área", y en grandes deformaciones se bifurca según si esa área es la previa o la posterior a la deformación. Una tabla lo resume.

Tensor	Cara donde actúa la fuerza	Área que divide	Simetría	Velocidad de deformación conjugada
Cauchy $\sigma$	deformada	deformada	simétrico	$d$ (pero como $J\sigma:d$ )
1er PK $P$	deformada	sin deformar	no simétrico	$\dot{F}$
2do PK $S$	retraída a la referencia	sin deformar	simétrico	$\dot{E}$
ingenieril $\varepsilon$ · $\sigma$	sin distinción	sin distinción	simétrico	solo en el límite infinitesimal

Las relaciones entre ellos son estas.

P = FS, \qquad \sigma = J^{-1} F S F^{T}, \qquad J = \det F

$J$ es la razón de volúmenes. El tensor 1er PK $P$ es no simétrico porque sus dos patas se apoyan en configuraciones distintas. Un índice apunta al estado deformado y el otro al no deformado. Por eso, guardar $P$ en un código de elementos finitos obliga a llevar las 9 componentes. $S$ pone ambas patas en la configuración de referencia y le bastan 6.

Que el par sea conjugado significa que la potencia coincide#

El "par conjugado" (conjugate pair) no es cuestión de gusto. Es la identidad de que la potencia interna por unidad de volumen de referencia debe dar el mismo valor.

\dot{W} = S : \dot{E} = P : \dot{F} = J\,\sigma : d

Aquí $d = \operatorname{sym}(\dot{F}F^{-1})$ es el tensor de velocidad de deformación (rate of deformation). Las tres expresiones escriben una misma magnitud física en tres configuraciones, así que sus valores deben coincidir exactamente.

En cambio, $\sigma : \dot{E}$ o $S : d$ no son magnitud física alguna. Las unidades cuadran y el cálculo corre, pero ese número no es una potencia. Si el residuo de elementos finitos se arma con esas combinaciones, la matriz de rigidez deja de ser el hessiano de un funcional de energía. Ensamblar ese residuo es el mismo trabajo que el método de rigidez directa de una cercha. Sin energía que minimizar, la iteración de Newton pierde la convergencia cuadrática.

Midiendo en Python la potencia de los tres pares en el mismo instante#

Se construyó una trayectoria de deformación que mezcla estiramiento, corte y rotación, y se midieron las tres combinaciones en $t=0.7$ . El material es Saint Venant-Kirchhoff, con $S = \lambda\,\mathrm{tr}(E)I + 2\mu E$ .

import math
 
I3 = [[1.0, 0, 0], [0, 1.0, 0], [0, 0, 1.0]]
 
def mul(A, B):
    return [[sum(A[i][k]*B[k][j] for k in range(3)) for j in range(3)] for i in range(3)]
 
def tr(A):
    return [[A[j][i] for j in range(3)] for i in range(3)]
 
def add(A, B, s=1.0):
    return [[A[i][j] + s*B[i][j] for j in range(3)] for i in range(3)]
 
def scale(A, s):
    return [[s*A[i][j] for j in range(3)] for i in range(3)]
 
def ddot(A, B):
    return sum(A[i][j]*B[i][j] for i in range(3) for j in range(3))
 
def trace(A):
    return A[0][0] + A[1][1] + A[2][2]
 
def det(A):
    return (A[0][0]*(A[1][1]*A[2][2] - A[1][2]*A[2][1])
          - A[0][1]*(A[1][0]*A[2][2] - A[1][2]*A[2][0])
          + A[0][2]*(A[1][0]*A[2][1] - A[1][1]*A[2][0]))
 
def inv(A):
    d = det(A)
    C = [[0.0]*3 for _ in range(3)]
    for i in range(3):
        for j in range(3):
            m = [[A[r][c] for c in range(3) if c != j] for r in range(3) if r != i]
            C[j][i] = ((-1)**(i+j))*(m[0][0]*m[1][1] - m[0][1]*m[1][0])/d
    return C
 
def sym(A):
    return scale(add(A, tr(A)), 0.5)
 
def green_lagrange(F):
    return scale(add(mul(tr(F), F), I3, -1.0), 0.5)
 
def linear_strain(F):
    return sym(add(F, I3, -1.0))
 
LAM, MU = 100.0, 60.0                      # constantes de Saint Venant-Kirchhoff
 
def pk2_stress(E):
    return add(scale(I3, LAM*trace(E)), E, 2*MU)
 
def cauchy_stress(F, S):
    return scale(mul(mul(F, S), tr(F)), 1.0/det(F))
 
def rot_z(th):
    c, s = math.cos(th), math.sin(th)
    return [[c, -s, 0.0], [s, c, 0.0], [0.0, 0.0, 1.0]]
 
print("--- 1. pure rotation, no stretch ---")
print(" theta   eps_xx(linear)   E_xx(Green-Lagrange)")
for deg in (0, 5, 10, 30, 60, 90):
    F = rot_z(math.radians(deg))
    print("%5.0f   %14.5f   %20.2e" % (deg, linear_strain(F)[0][0], green_lagrange(F)[0][0]))
 
def defo_path(t):
    """Aplica estiramiento y corte, luego una rotación rígida de 40 grados * t"""
    U = [[1 + 0.20*t, 0.15*t,     0.0],
         [0.15*t,     1 - 0.05*t, 0.0],
         [0.0,        0.0,        1 - 0.08*t]]
    return mul(rot_z(math.radians(40.0)*t), U)
 
def rate(f, t, h=1e-6):
    A, B = f(t + h), f(t - h)
    return [[(A[i][j] - B[i][j])/(2*h) for j in range(3)] for i in range(3)]
 
print()
print("--- 2. work rate at t=0.7, three pairings ---")
t = 0.7
F = defo_path(t)
Fd = rate(defo_path, t)
E = green_lagrange(F)
Ed = rate(lambda s: green_lagrange(defo_path(s)), t)
S = pk2_stress(E)
P = mul(F, S)
J = det(F)
sig = cauchy_stress(F, S)
d = sym(mul(Fd, inv(F)))                   # tensor de velocidad de deformación
 
print("  J = det F              = %.6f" % J)
print("  S : Edot   (2nd PK  x GL rate)   = %12.6f" % ddot(S, Ed))
print("  P : Fdot   (1st PK  x F rate)    = %12.6f" % ddot(P, Fd))
print("  J sigma: d (Cauchy  x stretching)= %12.6f" % (J*ddot(sig, d)))
print("  sigma : Edot   <- wrong pair     = %12.6f" % ddot(sig, Ed))
print("  S : d          <- wrong pair     = %12.6f" % ddot(S, d))

--- 1. pure rotation, no stretch ---
 theta   eps_xx(linear)   E_xx(Green-Lagrange)
    0          0.00000               0.00e+00
    5         -0.00381               0.00e+00
   10         -0.01519              -5.55e-17
   30         -0.13397               0.00e+00
   60         -0.50000               0.00e+00
   90         -1.00000               0.00e+00
 
--- 2. work rate at t=0.7, three pairings ---
  J = det F              = 1.028087
  S : Edot   (2nd PK  x GL rate)   =    10.522306
  P : Fdot   (1st PK  x F rate)    =    10.522306
  J sigma: d (Cauchy  x stretching)=    10.522306
  sigma : Edot   <- wrong pair     =     9.961790
  S : d          <- wrong pair     =     4.823514

5% y 54% — dos clases de error por equivocar el par#

Las tres combinaciones correctas coinciden hasta el sexto decimal. Dan $10.522306$ . Los números confirman que repartir la escritura en tres configuraciones no altera el valor.

Las dos combinaciones erróneas fallan de manera distinta. $\sigma : \dot{E}$ da $9.961790$ , un 5.3% por debajo. $\sigma$ y $S$ difieren por el mapeo $J^{-1}F(\cdot)F^{T}$ , y aquí ese mapeo es suave. Con $J = 1.028$ y un estiramiento de apenas un 20%, el error se detiene en ese orden.

$S : d$ da $4.823514$ . Un 54% por debajo. Este no es un problema de escala sino un error de otra especie. Ignora la relación $\dot{E} = F^{T} d F$ e introduce $d$ tal cual, de modo que también se mezclan las componentes de rotación. Cuanto mayor es el ángulo de rotación, mayor crece este error.

En la práctica el peligroso es el 5%. El 54% diverge en el primer paso de carga y se detecta de inmediato. El 5% sí converge, pero converge con la respuesta algo equivocada. Aumentar el número de elementos no lo elimina.

En la tabla de rotación pura la deformación ingenieril vale -0.00381 a 5°, -0.01519 a 10°, -0.13397 a 30°, -0.50000 a 60° y -1.00000 a 90°. Es $\cos\theta-1$ sin añadidos. La columna de Green-Lagrange es -5.55e-17 a 10° y cero en el resto. El -13% a 30° es un coseno, no un error de malla.

En $t=0.7$ los tres pares correctos coinciden en 10.522306 hasta seis decimales. $J=1.028087$ . $\sigma:\dot{E}=9.961790$ da $(10.522306-9.961790)/10.522306\approx 0.053$ , un 5.3% por debajo. $S:d=4.823514$ da $4.823514/10.522306\approx 0.458$ , así que falta cerca de un 54%. El 5% se esconde dentro de una respuesta ya convergida. El 54% revienta en el primer paso de carga.

El espesor lo cierra el volumen, no una ecuación#

Hasta aquí va la mecánica de medios continuos general. Las láminas suman un elemento más.

El elemento de lámina 3D lleva el estiramiento en el espesor como incógnita. En la formulación 3D-shell de Sussman y Bathe, la dirección del espesor necesita 3 incógnitas, mientras que la condición de tensión plana solo aporta 2 ecuaciones. Falta una.

Lo que llena la que falta es la condición de incompresibilidad.

J = \lambda_1 \lambda_2 \lambda_3 = 1 \quad\Longrightarrow\quad \lambda_3 = \frac{1}{\lambda_1 \lambda_2}

$\lambda_i$ son las razones de estiramiento principales. En materiales que casi conservan el volumen, como el caucho o la plasticidad metálica, el espesor deja de ser incógnita independiente y pasa a ser variable dependiente del estiramiento en el plano. Al estirar en el plano con $\lambda_1 = 1.20$ y $\lambda_2 = 1.05$ , el espesor queda en $0.794$ veces el original, un 20.6% más delgado. Es justo la relación que se usa para calcular el adelgazamiento en el conformado de chapa.

Hay un punto donde esta condición se encuentra con el atado MITC. Interpolar el estiramiento del espesor tal cual dentro del elemento genera bloqueo volumétrico (volumetric locking) en elementos delgados. Igual que el bloqueo por corte se resolvió con puntos de atado, el estiramiento del espesor se deja independiente solo en unos pocos puntos por elemento y el resto se ata por interpolación.

Truncar la rotación en primer orden hace crecer el director un 16%#

El segundo elemento propio de las láminas es la rotación. El elemento de lámina lleva consigo el vector normal a la superficie media, es decir, el director. Cuando la iteración de Newton entrega un incremento de rotación $\Delta\boldsymbol{\theta}$ , hay que girar el director esa cantidad. La matriz de rotación se construye con la fórmula de Rodrigues.

R(\boldsymbol{\theta}) = I + \frac{\sin\theta}{\theta}\,\Theta + \frac{1 - \cos\theta}{\theta^{2}}\,\Theta^{2}

$\Theta$ es la matriz antisimétrica de $\boldsymbol{\theta}$ y $\theta = |\boldsymbol{\theta}|$ . Es común el código que trunca a $R \approx I + \Theta$ alegando que el incremento es pequeño. Esa matriz no es ortogonal. Como $\det(I + \Theta) = 1 + \theta^{2} \neq 1$ , el director se alarga un poco en cada paso.

import math
 
def matvec(A, v):
    return [sum(A[i][k]*v[k] for k in range(3)) for i in range(3)]
 
def matmul(A, B):
    return [[sum(A[i][k]*B[k][j] for k in range(3)) for j in range(3)] for i in range(3)]
 
def skew(w):
    return [[0.0, -w[2], w[1]],
            [w[2], 0.0, -w[0]],
            [-w[1], w[0], 0.0]]
 
def rodrigues(w, order):
    """order = 1, 2 truncan la serie; 0 es la forma cerrada"""
    th = math.sqrt(sum(c*c for c in w))
    W = skew(w)
    W2 = matmul(W, W)
    if order == 1:
        a, b = 1.0, 0.0
    elif order == 2:
        a, b = 1.0, 0.5
    else:
        a = math.sin(th)/th
        b = (1.0 - math.cos(th))/(th*th)
    return [[(1.0 if i == j else 0.0) + a*W[i][j] + b*W2[i][j]
             for j in range(3)] for i in range(3)]
 
def spin_director(dth, steps, order):
    """Gira la normal de la superficie media alrededor del eje y, dth radianes por paso"""
    d = [0.0, 0.0, 1.0]
    for _ in range(steps):
        d = matvec(rodrigues([0.0, dth, 0.0], order), d)
    return d
 
print("--- 3. director after 30 increments of 0.10 rad (exact total 171.89 deg) ---")
print(" order        |d|      length err %   angle(deg)   angle err(deg)")
for order, name in ((1, "1st"), (2, "2nd"), (0, "closed")):
    d = spin_director(0.10, 30, order)
    n = math.sqrt(sum(c*c for c in d))
    ang = math.degrees(math.atan2(d[0], d[2]))
    if ang < 0:
        ang += 360.0
    print(" %-6s  %10.5f   %11.2f   %10.3f   %12.3f"
          % (name, n, 100*(n - 1.0), ang, ang - math.degrees(3.0)))
 
print()
print("--- 4. same total rotation, smaller increments (1st order) ---")
print(" steps   dtheta      |d|     length err %")
for steps in (30, 60, 150, 300, 3000):
    dth = 3.0/steps
    d = spin_director(dth, steps, 1)
    n = math.sqrt(sum(c*c for c in d))
    print(" %5d   %6.4f  %8.5f   %11.3f" % (steps, dth, n, 100*(n - 1.0)))
 
print()
print("--- 5. thickness closed by J = 1, not by a 3rd equation ---")
print(" lam1   lam2    lam3=1/(lam1 lam2)   thickness change %")
for l1, l2 in ((1.20, 1.05), (1.20, 1.00), (1.10, 1.10), (1.30, 0.95)):
    l3 = 1.0/(l1*l2)
    print(" %4.2f   %4.2f   %16.5f   %16.1f" % (l1, l2, l3, 100*(l3 - 1.0)))

--- 3. director after 30 increments of 0.10 rad (exact total 171.89 deg) ---
 order        |d|      length err %   angle(deg)   angle err(deg)
 1st        1.16097         16.10      171.318         -0.570
 2nd        1.00038          0.04      172.173          0.286
 closed     1.00000          0.00      171.887          0.000
 
--- 4. same total rotation, smaller increments (1st order) ---
 steps   dtheta      |d|     length err %
    30   0.1000   1.16097        16.097
    60   0.0500   1.07778         7.778
   150   0.0200   1.03045         3.045
   300   0.0100   1.01511         1.511
  3000   0.0010   1.00150         0.150
 
--- 5. thickness closed by J = 1, not by a 3rd equation ---
 lam1   lam2    lam3=1/(lam1 lam2)   thickness change %
 1.20   1.05            0.79365              -20.6
 1.20   1.00            0.83333              -16.7
 1.10   1.10            0.82645              -17.4
 1.30   0.95            0.80972              -19.0

Con d.theta en 0.10 la flecha roja de primer orden sale del círculo unitario. A la mitad, la deriva también se reduce a la mitad.

d.theta0.100

Leave d.theta at 0.10 rad and let it march: after 0 increments the red 1st-order director is 0.0 % too long and has climbed off the dashed unit circle, while the green closed-form arrow sits on it. Drag d.theta down — the drift shrinks in proportion, so halving the load step only halves the error. The yellow 2nd-order curve (0.00 %) shows what one more term buys.

A los 30 pasos, 171.89° en total, la longitud roja crece hacia el 1.16097 impreso, un error del 16.10%. La forma cerrada verde permanece sobre el círculo con error angular nulo.

El truncamiento de primer orden alargó el director un 16.1% en apenas 30 pasos. Con un término más, el segundo orden baja a 0.04%. La longitud ganó 400 veces en precisión, pero el error de ángulo resulta mayor justamente en segundo orden. Son dos errores independientes.

La tabla 4 importa más. Reducir el incremento a la décima parte reduce el error solo a la décima parte. Partir la carga en pasos más finos no elimina este problema. Hay que usar la forma cerrada, normalizar el director en cada paso, o llevar la rotación como cuaternión.

Con $\lambda_1=1.20$ y $\lambda_2=1.05$ , $\lambda_3=0.79365$ , un 20.6% más delgado. Con $\lambda_1=\lambda_2=1.10$ el espesor queda en 0.82645, un 17.4% menos. Partir la misma rotación en 3000 pasos de 0.001 rad deja todavía un error de longitud del 0.150%. Cien veces más fino, el error solo baja cerca de un factor cien. El segundo orden en los mismos 30 pasos tiene un error de longitud del 0.04% pero un error angular de +0.286°, peor que el 0.000° de la forma cerrada.

Si a 30° la tensión no es cero, $F$ no se elevó al cuadrado#

Si $B$ se construye como $\partial E / \partial u$ , la tensión que multiplica es $S$ . Si es $\partial \varepsilon / \partial u$ , es $\sigma$ . La definición de $B$ decide qué tensión es, el mismo punto que al cambiar de coordenadas. Integrar $\int_{V_0}$ sin un factor $J$ significa $S$ . Multiplicar por $J$ es retraer $\sigma$ a la referencia. Si $\sigma = J^{-1}FSF^{T}$ aparece justo antes de von Mises, el cálculo corría con $S$ . Saltar ese mapa y graficar $S$ directo no se nota hasta pasar un 20% de estiramiento. En esta trayectoria, con $J = 1.028$ , $\sigma:\dot{E}$ ya está un 5.3% por debajo.

Si esos sitios no se leen, se gira un solo elemento y nada más. A 30° la deformación ingenieril es -0.134, cerca de un 13%, y una tensión distinta de cero significa que $F^{T}F$ no se formó. Green-Lagrange es cero en ese ángulo. Un Rodrigues de primer orden alarga el director un 16.1% en 30 pasos de 0.10 rad. Partir la carga en 3000 pasos deja todavía un error de longitud del 0.15%.