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.

Conviene manipularlo directamente en la simulación de abajo.

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.

Al pulsar rigid rotation only y dejar correr la rotación, la forma del elemento no cambia, pero la barra roja de $\varepsilon$ oscila con fuerza de un lado a otro. La barra verde de $E$ permanece pegada al cero. Al subir lam, recién ahí se mueve $E$ , y su valor no cambia por más que se varíe el ángulo de rotación.

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. Es el mismo terreno que la simetría de Galerkin. 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.

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.

Conviene manipularlo directamente en la simulación de abajo.

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.

Con d.theta en 0.10 y la simulación en marcha, la flecha roja truncada a primer orden escapa en espiral fuera del círculo unitario punteado. La forma cerrada verde se mantiene exactamente sobre el círculo. Al reducir d.theta a la mitad, la deriva también se reduce a la mitad: el error es de primer orden en el tamaño del incremento.

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

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.

Cómo distinguir en el código si es $S$ o $\sigma$ #

Al abrir el código de grandes deformaciones de otra persona, la identidad de la variable de tensión no se deduce del nombre. Tres lugares la delatan.

Se mira el sitio donde la tensión entra en la matriz de rigidez. Si la matriz $B$ está construida como $\partial E / \partial u$ , la tensión que la multiplica es $S$ . Si $B$ es $\partial \varepsilon / \partial u$ , es $\sigma$ . Como se vio en la transformación de la matriz constitutiva, la definición de $B$ decide la identidad de la tensión.

Se mira el jacobiano de la integración. Si se integra sobre el volumen de referencia como $\int_{V_0} \cdots \, dV_0$ sin multiplicar por $J$ , es $S$ . Si se está multiplicando por $J$ , entonces se está retrayendo $\sigma$ a la configuración de referencia.

Se mira la rutina de salida. Si en el posproceso aparece la transformación $\sigma = J^{-1}FSF^{T}$ justo antes de calcular von Mises, es señal de que por dentro todo corría con $S$ . Existen códigos que se saltan esa transformación y grafican las componentes de $S$ metidas directamente en von Mises. En deformaciones pequeñas no se nota; pasado un 20% de estiramiento, las curvas se separan.

Si ninguna de las tres se puede confirmar, queda la prueba. Se somete un solo elemento a rotación rígida pura y se verifica que la tensión se mantenga en cero. Si a 30° aparece un 13%, es que en algún lado no se elevó $F$ al cuadrado.