Un mismo valor integral leído como 0.1406 y como 56.91 — función de tensión de torsión y flujo laminar en ductos

Basta resolver una vez el problema con laplaciano igual a -1 sobre la sección: ese valor integral es a la vez la constante de torsión y el f·Re del ducto. El solver estructural y el de flujo ensamblan la misma matriz dos veces.

Torcer una barra cuadrada y hacer correr agua por un ducto cuadrado#

¿Cuánto torque cuesta torcer 0.01 rad por metro una barra de acero de sección cuadrada? ¿Y cuánto vale el factor de fricción si por un ducto igual circula agua? Las dos preguntas se enseñan en facultades distintas y con libros distintos. Pero la respuesta sale de un único valor integral.

Este artículo calcula ese valor de forma directa. Se resuelve una vez la ecuación de Poisson sobre la sección con elementos triangulares lineales y esa misma solución se lee dos veces. Una vez como constante de torsión $J$ , otra como grupo de fricción laminar $f\cdot Re$ . En sección cuadrada los valores que deben aparecer son 0.1406 y 56.91. Ambos están en los manuales.

Cómo un problema tridimensional se reduce a un escalar sobre la sección#

Primero la torsión. Prandtl no resolvió las componentes de tensión de manera directa. En su lugar introdujo una función de tensión (stress function, un campo escalar cuya derivada entrega las tensiones) $\Phi$ .

\tau_{zx} = \frac{\partial \Phi}{\partial y}, \qquad \tau_{zy} = -\frac{\partial \Phi}{\partial x}

$\tau_{zx}$ y $\tau_{zy}$ son las dos componentes de la tensión cortante sobre la sección. Así planteado, el equilibrio queda satisfecho por construcción. Solo resta una condición de compatibilidad, y esa se convierte en la ecuación de Poisson sobre la sección $A$ .

\nabla^2 \Phi + 2G\theta = 0, \qquad \Phi = 0 \ \text{ on } \ \partial A

$G$ es el módulo de corte y $\theta$ el ángulo de torsión por unidad de longitud. La cara lateral es libre y no soporta corte, así que $\Phi$ es constante en el contorno. En sección maciza esa constante se puede fijar en 0. El torque se recupera con una integral de sección.

T = 2\int_A \Phi \, dA

Ahora el flujo. Al desarrollarse por completo en un conducto de sección constante, de la velocidad solo sobrevive la componente axial $u(x,y)$ . Como $u$ no cambia a lo largo del eje, el término convectivo desaparece entero. Lo que resta de Navier–Stokes es lineal.

\mu \nabla^2 u = \frac{dp}{dz}, \qquad u = 0 \ \text{ on the wall}

$\mu$ es la viscosidad dinámica y $dp/dz$ el gradiente axial de presión, constante sobre la sección. El caudal también es una integral de la misma forma.

Q = \int_A u \, dA

Las dos expresiones solo difieren en los símbolos. Conviene manipularlo en la simulación de abajo.

aspect ratio w/h1.0

J/(w h³) 0.00000f·Re 0.000sweeps 0

Drag the aspect ratio: the solve restarts and both cards refill from the one field. Watch that the two numbers move in opposite directions — J/(w h³) climbs toward 1/3 while f·Re climbs toward 96 — and that neither card ever needs a second matrix.

Al mover el deslizador de relación de aspecto, la relajación corre otra vez desde cero y las dos tarjetas de la derecha se llenan desde el mismo campo. Los botones no cambian el cálculo. Solo cambian la etiqueta.

Tabla de correspondencia, símbolo por símbolo#

Torsión	Flujo laminar en ducto	En común
Función de tensión $\Phi$	Velocidad axial $u$	Escalar incógnito
$2G\theta$	$(-dp/dz)/\mu$	Lado derecho constante
Superficie libre $\Phi = 0$	No deslizamiento $u = 0$	Frontera de Dirichlet
Tensión cortante $\lvert\nabla\Phi\rvert$	Corte en pared $\mu\lvert\nabla u\rvert$	Pendiente en la frontera
Torque $T = 2\int \Phi\, dA$	Caudal $Q = \int u\, dA$	Integral de sección
Constante de torsión $J$	$f\cdot Re$	Constante fijada por la forma

Conviene normalizar una sola vez. Se resuelve el problema $\nabla^2 w = -1$ con $w = 0$ en el contorno y se define $I = \int_A w \, dA$ . Entonces las dos constantes caen así.

J = 4I, \qquad f\cdot Re = \frac{2 D_h^2}{I/A}

La primera expresión sale de llevar $\Phi = 2G\theta\, w$ a $T = G\theta J$ . La segunda resulta de multiplicar el factor de fricción de Darcy $f = 2 D_h(-dp/dz)/(\rho \bar u^2)$ por $Re = \rho \bar u D_h/\mu$ y cancelar la velocidad media $\bar u$ . $D_h = 4A/P$ es el diámetro hidráulico y $P$ el perímetro mojado.

Meter la sección circular sirve de verificación. Con radio $R$ se tiene $w = (R^2 - r^2)/4$ , $I/A = R^2/8$ y $D_h = 2R$ . Al sustituir aparece $f\cdot Re = 64$ tal cual. De dónde sale el exponente 4 en el tubo circular ya se trató por separado.

El 3×3 que produce un solo triángulo lineal#

La forma débil es idéntica en ambos lados. Al multiplicar por la función de prueba e integrar por partes, en la matriz de rigidez solo queda el producto interno de los gradientes de las funciones de forma. En un triángulo lineal ese gradiente es constante dentro del elemento. No hacen falta puntos de integración y la matriz sale en forma cerrada.

K^{(e)}_{mn} = \frac{\beta_m \beta_n + \delta_m \delta_n}{4 A_e}, \qquad F^{(e)}_m = \frac{A_e}{3}

Para los nodos $i,j,k$ se tiene $\beta_i = Y_j - Y_k$ y $\delta_i = X_k - X_j$ ; el resto se obtiene rotando los índices. $A_e$ es el área del triángulo. El término de carga se reparte por igual en un tercio del área a cada nodo porque el lado derecho es constante. El camino que llega a la misma matriz por residuos ponderados de Galerkin quedó ordenado antes.

Python — un solo CG para dos constantes#

El código de abajo corta la sección rectangular con una malla estructurada y pone dos triángulos por celda. Resuelve una vez con CG precondicionado por la diagonal y lee la solución dos veces. La solución en serie sirve de verificación.

import math
 
 
def tri_stiffness(p0, p1, p2):
    """Linear triangle: K = (beta_i beta_j + delta_i delta_j) / (4A)."""
    (x0, y0), (x1, y1), (x2, y2) = p0, p1, p2
    a2 = x0 * (y1 - y2) + x1 * (y2 - y0) + x2 * (y0 - y1)
    area = 0.5 * a2
    beta = (y1 - y2, y2 - y0, y0 - y1)
    delta = (x2 - x1, x0 - x2, x1 - x0)
    k = [[(beta[r] * beta[c] + delta[r] * delta[c]) / (2.0 * a2)
          for c in range(3)] for r in range(3)]
    return k, area
 
 
def build_mesh(w, h, nx, ny):
    nodes, idx = [], {}
    for j in range(ny + 1):
        for i in range(nx + 1):
            idx[(i, j)] = len(nodes)
            nodes.append((w * i / nx, h * j / ny))
    tris = []
    for j in range(ny):
        for i in range(nx):
            a, b = idx[(i, j)], idx[(i + 1, j)]
            c, d = idx[(i + 1, j + 1)], idx[(i, j + 1)]
            tris.append((a, b, c))
            tris.append((a, c, d))
    fixed = set()
    for j in range(ny + 1):
        for i in range(nx + 1):
            if i in (0, nx) or j in (0, ny):
                fixed.add(idx[(i, j)])
    return nodes, tris, fixed
 
 
def assemble_poisson(nodes, tris, fixed):
    """-lap(u) = 1 with u = 0 on 'fixed'.  Returns CSR-ish rows and rhs."""
    n = len(nodes)
    rows = [dict() for _ in range(n)]
    rhs = [0.0] * n
    for (a, b, c) in tris:
        k, area = tri_stiffness(nodes[a], nodes[b], nodes[c])
        ids = (a, b, c)
        for r in range(3):
            if ids[r] in fixed:
                continue
            rhs[ids[r]] += area / 3.0
            for c2 in range(3):
                if ids[c2] in fixed:
                    continue
                rows[ids[r]][ids[c2]] = rows[ids[r]].get(ids[c2], 0.0) + k[r][c2]
    for f in fixed:
        rows[f] = {f: 1.0}
        rhs[f] = 0.0
    return rows, rhs
 
 
def cg_solve(rows, rhs, tol=1e-12, itmax=20000):
    n = len(rhs)
    x = [0.0] * n
    r = rhs[:]
    z = [r[i] / rows[i][i] for i in range(n)]
    p = z[:]
    rz = sum(r[i] * z[i] for i in range(n))
    r0 = math.sqrt(sum(v * v for v in r))
    for it in range(itmax):
        ap = [0.0] * n
        for i in range(n):
            s = 0.0
            for j, v in rows[i].items():
                s += v * p[j]
            ap[i] = s
        alpha = rz / sum(p[i] * ap[i] for i in range(n))
        for i in range(n):
            x[i] += alpha * p[i]
            r[i] -= alpha * ap[i]
        rn = math.sqrt(sum(v * v for v in r))
        if rn <= tol * r0:
            return x, it + 1
        z = [r[i] / rows[i][i] for i in range(n)]
        rz2 = sum(r[i] * z[i] for i in range(n))
        beta = rz2 / rz
        rz = rz2
        p = [z[i] + beta * p[i] for i in range(n)]
    return x, itmax
 
 
def section_integral(nodes, tris, u):
    tot = 0.0
    for (a, b, c) in tris:
        _, area = tri_stiffness(nodes[a], nodes[b], nodes[c])
        tot += area * (u[a] + u[b] + u[c]) / 3.0
    return tot
 
 
def series_rect(w, h, nterm=60):
    """Exact integral of the Prandtl/duct solution over a w x h rectangle."""
    s = h if h < w else w
    lg = w if h < w else h
    acc = 0.0
    for m in range(1, 2 * nterm, 2):
        acc += math.tanh(m * math.pi * lg / (2.0 * s)) / m ** 5
    j = (1.0 / 3.0) * lg * s ** 3 * (1.0 - (192.0 / math.pi ** 5) * (s / lg) * acc)
    return j / 4.0          # integral of u  ==  J / 4
 
 
def solve_section(w, h, nx, ny):
    nodes, tris, fixed = build_mesh(w, h, nx, ny)
    rows, rhs = assemble_poisson(nodes, tris, fixed)
    u, its = cg_solve(rows, rhs)
    iu = section_integral(nodes, tris, u)
    area, perim = w * h, 2.0 * (w + h)
    dh = 4.0 * area / perim
    return dict(int_u=iu, jtor=4.0 * iu, umean=iu / area,
                fre=2.0 * dh * dh / (iu / area), umax=max(u), its=its,
                ndof=len(nodes))
 
 
if __name__ == '__main__':
    ex_i = series_rect(1.0, 1.0)
    print("[A] square bar, one Poisson solve -> two constants  (exact J/a^4 = %.6f,"
          " f*Re = %.4f)" % (4 * ex_i, 2.0 / ex_i))
    print("  mesh    nodes    integral u      J/a^4    err(%)      f*Re   err(%)  CG")
    for n in (8, 16, 32, 64):
        r = solve_section(1.0, 1.0, n, n)
        print("%3dx%-3d %7d  %.8f  %.6f  %7.3f  %8.4f  %7.3f %4d"
              % (n, n, r['ndof'], r['int_u'], r['jtor'],
                 100 * (r['jtor'] / (4 * ex_i) - 1), r['fre'],
                 100 * (r['fre'] / (2.0 / ex_i) - 1), r['its']))
 
    print()
    print("[B] aspect-ratio sweep, 64x64 mesh   (w x h = AR x 1)")
    print("  AR   J/(w h^3)   beta2(ref)     f*Re   fRe(exact)   u_max/u_mean")
    REF_B2 = {1: 0.1406, 2: 0.2290, 4: 0.2810, 8: 0.3070}
    for ar in (1, 2, 4, 8):
        w, h = float(ar), 1.0
        r = solve_section(w, h, 64, 64)
        ex = series_rect(w, h)
        dh = 4.0 * w * h / (2.0 * (w + h))
        print("%4d   %.5f      %.4f    %8.4f   %8.4f   %8.4f"
              % (ar, r['jtor'] / (w * h ** 3), REF_B2[ar], r['fre'],
                 2 * dh * dh / (ex / (w * h)), r['umax'] / r['umean']))
 
    print()
    print("[C] the same integral read twice (square section, 64x64)")
    r = solve_section(1.0, 1.0, 64, 64)
    print("  dimensionless integral of u over A = %.8f a^4" % r['int_u'])
 
    G, THETA, SIDE = 80e9, 0.01, 0.05          # steel bar, 50 mm square
    j = r['jtor'] * SIDE ** 4
    print("  steel bar  a=50 mm, G=80 GPa, twist=0.01 rad/m")
    print("     J = 4*int*a^4 = %.4e m^4    T = G*theta*J = %.1f N.m" % (j, G * THETA * j))
 
    MU, RHO, DPDX, HALF = 1.0e-3, 1000.0, 200.0, 0.005   # water, 5 mm square duct
    q = DPDX / MU * r['int_u'] * HALF ** 4
    area = HALF ** 2
    ubar = q / area
    dh = HALF
    re = RHO * ubar * dh / MU
    print("  water duct a=5 mm, dp/dx=200 Pa/m, mu=1e-3 Pa.s")
    print("     Q = (G_p/mu)*int*a^4 = %.3e m^3/s   u_mean = %.4f m/s   Re = %.0f"
          % (q, ubar, re))
    print("     f = (f*Re)/Re = %.4f      f*Re = %.4f  (handbook 56.91)"
          % (r['fre'] / re, r['fre']))
    print("  J/a^4 (bar) = %.8f   vs   4*mu*Q/(G_p*a^4) (duct) = %.8f"
          % (j / SIDE ** 4, 4 * MU * q / DPDX / HALF ** 4))

[A] square bar, one Poisson solve -> two constants  (exact J/a^4 = 0.140577, f*Re = 56.9083)
  mesh    nodes    integral u      J/a^4    err(%)      f*Re   err(%)  CG
  8x8        81  0.03342303  0.133692   -4.898   59.8390    5.150    9
 16x16      289  0.03470275  0.138811   -1.256   57.6323    1.272   32
 32x32     1089  0.03503302  0.140132   -0.317   57.0890    0.318   70
 64x64     4225  0.03511638  0.140466   -0.079   56.9535    0.079  142
 
[B] aspect-ratio sweep, 64x64 mesh   (w x h = AR x 1)
  AR   J/(w h^3)   beta2(ref)     f*Re   fRe(exact)   u_max/u_mean
   1   0.14047      0.1406     56.9535    56.9083     2.0975
   2   0.22848      0.2290     62.2469    62.1922     1.9932
   4   0.28047      0.2810     73.0194    72.9311     1.7758
   8   0.30647      0.3070     82.4998    82.3386     1.6315
 
[C] the same integral read twice (square section, 64x64)
  dimensionless integral of u over A = 0.03511638 a^4
  steel bar  a=50 mm, G=80 GPa, twist=0.01 rad/m
     J = 4*int*a^4 = 8.7791e-07 m^4    T = G*theta*J = 702.3 N.m
  water duct a=5 mm, dp/dx=200 Pa/m, mu=1e-3 Pa.s
     Q = (G_p/mu)*int*a^4 = 4.390e-06 m^3/s   u_mean = 0.1756 m/s   Re = 878
     f = (f*Re)/Re = 0.0649      f*Re = 56.9535  (handbook 56.91)
  J/a^4 (bar) = 0.14046553   vs   4*mu*Q/(G_p*a^4) (duct) = 0.14046553

La última línea de [C] es el punto del artículo. El $J/a^4$ de la barra de acero de 50 mm y el $4\mu Q/(G_p a^4)$ del ducto de agua de 5 mm son el mismo número hasta el octavo decimal. Uno entrega 702.3 N·m; el otro, 4.39 mL por segundo.

En qué proporción se mueve el error al reducir el elemento a la mitad#

El error de [A] va 4.898 → 1.256 → 0.317 → 0.079 %. Cada vez que la malla se reduce a la mitad, baja alrededor de 4 veces. La solución con triángulos lineales es $O(h^2)$ y su integral sigue el mismo orden.

El signo resulta más interesante. Las cuatro mallas ven $J$ por debajo y $f\cdot Re$ por encima. No es casualidad. La solución de desplazamientos por elementos finitos siempre es más rígida que la exacta. Con la rigidez sobreestimada, el mismo torque tuerce menos y el valor integral baja. Traducido al flujo, el caudal se subestima. Con menos caudal, el factor de fricción sube. Es un caso poco común en el que el mismo sesgo se lee del lado seguro en ambos lados.

Las iteraciones de CG crecieron 9 → 32 → 70 → 142. Van casi en proporción al número de nodos por lado de la malla. Es el comportamiento típico de Poisson, con número de condición que crece como $h^{-2}$ , y el punto en el que una sección más grande empieza a pedir multigrid.

Al forzar la relación de aspecto, los caminos se separan en 0.333 y 96#

En [B], los $J/(wh^3)$ para relaciones de aspecto 1, 2, 4 y 8 son 0.1405, 0.2285, 0.2805 y 0.3065. La tabla de $\beta_2$ de los textos de mecánica de materiales da 0.1406, 0.229, 0.281 y 0.307, así que coinciden hasta el tercer decimal. En las mismas filas, $f\cdot Re$ vale 56.95, 62.25, 73.02 y 82.50. La solución en serie da 56.91, 62.19, 72.93 y 82.34.

Las dos constantes crecen en paralelo, pero sus límites son distintos. Cuando la sección se adelgaza, $J/(wh^3)$ tiende a $1/3$ y $f\cdot Re$ tiende al 96 de las placas paralelas. Con relación de aspecto 8 ya se llegó a 0.3065 y 82.5.

La última columna es la razón entre velocidad máxima y velocidad media. En la sección cuadrada salió 2.0975 y el valor de manual es 2.096. Cuanto más delgada, más baja hacia el 1.5 de las placas paralelas. Con relación de aspecto 8 es 1.63. Esta columna no tiene contraparte del lado de la torsión. En estructuras interesa la tensión máxima; en flujo, la velocidad máxima.

La película de jabón señala dónde está la tensión máxima#

Prandtl también dejó una manera de leer esta ecuación por experimento, sin calcularla. Si se cubre con una película de jabón un agujero con la forma de la sección y se sopla suavemente, la deflexión de la película es $\Phi$ , su pendiente es la tensión cortante y el volumen que desplaza es el torque. La membrana bajo presión uniforme obedece la misma ecuación de Poisson.

aspect ratio w/h2.0

peak on the longc1 0.0000τmax/τmean 0.000

Stretch the section and watch the red marker: it stays at the middle of the long side and never moves to a corner, where the film is flat. Switch the two buttons — the bars are not redrawn, only relabelled.

Conviene seguir el punto rojo mientras se alarga la relación de aspecto. La pendiente máxima se sienta siempre en el centro del lado largo y se apaga a 0 cerca de las esquinas. En la esquina la película queda sujeta por dos lados a la vez y está casi plana.

En la práctica esta sola línea rinde bastante. En una barra cuadrada torsionada, la grieta arranca en el centro del lado largo, no en la esquina. En un ducto, el corte en pared es máximo en ese mismo centro y casi nulo en el rincón. Que se asienten sedimentos en las esquinas de un ducto cuadrado, y que los productos de corrosión del rincón no se laven bien, son la misma figura. Cuanto más delgada la sección, más plano queda el corte en pared y más se acerca a 1 la razón entre máximo y promedio.

Dónde se rompe esta correspondencia#

Se rompe primero en secciones huecas. Con más de un contorno, $\Phi$ toma una constante distinta en cada uno y aparecen condiciones extra que fijan esas constantes. Del lado del flujo no existe nada parecido. Es apenas un problema de Dirichlet con una pared más.

Del lado del flujo caen primero las hipótesis. En la región de entrada $u$ cambia a lo largo del eje y el término convectivo revive; pasado $Re$ de 2000, la frase misma de que $f\cdot Re$ es constante deja de sostenerse. Si la viscosidad queda arrastrada por la temperatura, o si entran superficie libre y flotabilidad, el lado derecho ya no es constante sobre la sección.

Del lado de la torsión el papel equivalente lo cumplen la plasticidad y la restricción al alabeo. Una sección abierta delgada con los extremos sujetos queda fuera de la hipótesis de St. Venant.

Las condiciones que quedan caben en dos líneas. Que el lado derecho sea constante sobre la sección. Que toda la frontera sea de Dirichlet. Mientras esas dos sigan vivas, los dos problemas son el mismo problema.

Se estaba ensamblando la misma matriz dos veces#

El módulo de torsión de un código estructural y el módulo de flujo desarrollado de un código de flujo son la misma rutina de ensamblaje escrita dos veces. Solo cambia una constante en el lado derecho y la etiqueta que se pega al valor integral una vez resuelto.

Por eso la verificación también se despacha de una sola vez. Si hay un solver de torsión recién escrito, basta meterle la sección cuadrada y sacar también el $f\cdot Re$ . Si aparece 56.91, el 0.1406 del lado estructural también está bien. Es el mismo número, no hay margen de error.