회전만 30° 시켰는데 변형률이 -13%로 나왔다 — 대변형 쉘의 응력·변형률 짝짓기

대변형에서는 응력과 변형률을 아무거나 곱하면 안 된다. $S:\dot{E}$, $P:\dot{F}$, $J\sigma:d$만 같은 값을 낸다.

회전만 시켰는데 변형률 게이지가 -13%를 가리켰다#

쉘 요소 하나를 평면 안에서 30° 돌렸다. 늘이지도 않았고 비틀지도 않았다. 강체 회전만 했다.

그런데 공학 변형률 $\varepsilon_{xx} = \partial u_x / \partial x$ 를 계산하면 $-0.134$ 가 나온다. 13.4% 압축이다. 요소는 아무 데도 변형되지 않았는데 게이지는 압축을 읽는다.

이 값은 실수도 아니고 이산화 오차도 아니다. $\cos 30° - 1 = -0.134$ , 정확히 그 값이다. 미소변형 변형률의 정의 자체가 회전을 변형으로 오독한다.

이 글은 그 오독을 어디서 끊는지, 그리고 끊고 나면 응력을 어떤 것으로 바꿔 들어야 하는지를 다룬다. 얻어 갈 것은 세 가지다. 회전에 면역인 변형률 척도, 응력·변형률 짝을 틀렸을 때 에너지가 얼마나 어긋나는지의 실측값, 그리고 쉘 두께 방향 미지수를 닫는 방정식의 정체다.

아래 시뮬레이션에서 직접 조작해보자.

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.

rigid rotation only을 누르고 회전을 흘려보자. 요소 모양은 그대로인데 빨간 $\varepsilon$ 막대만 좌우로 크게 흔들린다. 초록 $E$ 막대는 0에 붙어 있다. lam을 올리면 그제야 $E$ 가 움직이고, 그 값은 회전각을 아무리 바꿔도 변하지 않는다.

회전을 걸러내는 것은 $F$ 가 아니라 $F^{T}F$ 다#

변형은 변형 구배(deformation gradient, 기준 배치에서 현재 배치로의 국소 사상)로 시작한다.

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

$x$ 는 변형 후 좌표, $X$ 는 변형 전 좌표다. 극분해(polar decomposition)를 쓰면 $F = RU$ 로 갈라진다. $R$ 은 회전, $U$ 는 순수 늘임이다. 문제는 $F$ 자체에 $R$ 이 살아 있다는 것이다.

$\varepsilon = \tfrac{1}{2}(F + F^{T}) - I$ 는 $F$ 를 그대로 쓴다. 그래서 $R$ 이 새어 들어온다. 회전만 있을 때 $\varepsilon_{xx} = \cos\theta - 1$ 이 되는 이유가 이것이다.

그런데 $F$ 를 제곱하면 회전이 사라진다.

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

$R^{T}R = I$ 이므로 $C$ 에는 $U$ 만 남는다. 이 $C$ 가 우 Cauchy-Green 텐서다. 여기서 항등을 빼고 반으로 나눈 것이 Green-Lagrange 변형률이다.

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

강체 회전에서 $F^{T}F = I$ 이므로 $E$ 는 정확히 0이다. 위 시뮬레이션의 초록 막대가 움직이지 않는 이유가 이 두 줄이다. 좌표계를 바꿔도 물리가 변하지 않아야 한다는 요구를 구성 텐서 좌표변환에서는 기저 변환으로 풀었다면, 여기서는 변형률 척도의 정의 자체로 푼다.

응력 텐서는 어느 면적으로 나눈 값인가#

변형률을 기준 배치로 끌고 왔으면 응력도 같이 끌고 와야 한다. 응력은 "힘 나누기 면적"인데, 대변형에서는 그 면적이 변형 전인지 후인지가 갈린다. 표 한 장으로 정리된다.

텐서	힘이 작용하는 면	나누는 면적	대칭성	짝이 되는 변형률 속도
Cauchy $\sigma$	변형 후	변형 후	대칭	$d$ (단, $J\sigma:d$ )
1st PK $P$	변형 후	변형 전	비대칭	$\dot{F}$
2nd PK $S$	기준 배치로 당김	변형 전	대칭	$\dot{E}$
공학 $\varepsilon$ · $\sigma$	구분 없음	구분 없음	대칭	미소변형 극한에서만

서로의 관계는 다음과 같다.

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

$J$ 는 부피비다. 1st PK 텐서 $P$ 가 비대칭인 이유는 두 다리가 서로 다른 배치를 딛고 있기 때문이다. 한쪽 첨자는 변형 후, 다른 쪽은 변형 전을 가리킨다. 그래서 유한요소 코드에서 $P$ 를 저장하려면 9개 성분을 다 들고 있어야 한다. $S$ 는 양쪽 다리를 기준 배치에 놓아 6개면 된다.

짝이 맞는다는 것은 일률이 같다는 뜻이다#

"짝(conjugate pair)"은 취향 문제가 아니다. 단위 기준 부피당 내부 일률이 같은 값을 내야 한다는 등식이다.

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

여기서 $d = \operatorname{sym}(\dot{F}F^{-1})$ 는 변형률 속도(rate of deformation) 텐서다. 세 표현은 같은 물리량을 세 배치에서 적은 것이므로 값이 정확히 일치해야 한다.

반대로 $\sigma : \dot{E}$ 나 $S : d$ 는 아무 물리량도 아니다. 단위는 맞고 계산도 되지만 그 수는 일률이 아니다. 유한요소 잔차를 이런 조합으로 세우면 강성행렬이 에너지 범함수의 헤시안이 되지 못한다. 갈러킨 대칭성 이야기와 같은 자리다. 최소화할 에너지가 없으면 뉴턴 반복은 이차 수렴을 잃는다.

Python으로 세 짝의 일률을 같은 시각에 재봤다#

늘임과 전단과 회전을 섞은 변형 경로를 하나 만들고, $t=0.7$ 에서 세 조합의 값을 재봤다. 재료는 Saint Venant-Kirchhoff, $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                      # 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):
    """늘임과 전단을 준 뒤 40도 * 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)))                   # 변형률 속도 텐서
 
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%와 54% — 짝을 틀리면 생기는 두 종류의 오차#

정상인 세 조합은 소수점 여섯 자리까지 같다. $10.522306$ 이다. 배치를 셋으로 나눠 적었을 뿐 같은 값이라는 것이 숫자로 확인된다.

틀린 조합 둘은 서로 다르게 틀린다. $\sigma : \dot{E}$ 는 $9.961790$ 으로 5.3% 낮다. $\sigma$ 와 $S$ 가 $J^{-1}F(\cdot)F^{T}$ 만큼 다른데 이 변형이 작기 때문이다. $J = 1.028$ 이고 늘임도 20% 남짓이라 오차도 그 정도에서 멈춘다.

$S : d$ 는 $4.823514$ 다. 54% 낮다. 이쪽은 스케일 문제가 아니라 종류가 다른 오류다. $\dot{E} = F^{T} d F$ 라는 관계를 무시하고 $d$ 를 그대로 넣었기 때문에 회전 성분까지 섞여 들어간다. 회전각이 커질수록 이 오차는 커진다.

실무에서 위험한 쪽은 5%다. 54%는 첫 하중 스텝에서 발산해 바로 잡힌다. 5%는 수렴은 하는데 답이 조금 틀린 채로 수렴한다. 요소 개수를 늘려도 사라지지 않는다.

두께 방향은 방정식이 아니라 부피가 닫는다#

여기까지가 일반 연속체 이야기다. 쉘에는 항목이 하나 더 붙는다.

3D 쉘 요소는 두께 방향 늘임을 미지수로 들고 있다. Sussman과 Bathe의 3D-shell 정식화에서는 두께 방향에 미지수가 3개 필요한데, 평면응력 조건에서 나오는 식은 2개뿐이다. 하나가 모자란다.

모자란 하나를 채우는 것이 비압축 조건이다.

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

$\lambda_i$ 는 주늘임비다. 고무나 금속 소성처럼 부피가 거의 보존되는 재료에서, 두께는 독립 미지수가 아니라 면내 늘임의 종속변수가 된다. 면내로 $\lambda_1 = 1.20$ , $\lambda_2 = 1.05$ 만큼 늘이면 두께는 $0.794$ 배, 20.6% 얇아진다. 판재 성형에서 두께 감소를 계산할 때 쓰는 바로 그 관계다.

이 조건이 MITC 타잉과 만나는 지점이 있다. 두께 방향 늘임을 요소 내부에서 그대로 보간하면 얇은 요소에서 체적 잠김(volumetric locking)이 생긴다. 전단 잠김을 타잉점으로 풀었듯, 두께 늘임도 요소당 몇 개의 점에서만 독립으로 두고 나머지는 보간으로 묶는다.

회전을 1차에서 끊으면 director가 16% 자란다#

쉘의 두 번째 항목은 회전이다. 쉘 요소는 중립면 법선 벡터, 즉 director를 들고 다닌다. 뉴턴 반복이 회전 증분 $\Delta\boldsymbol{\theta}$ 를 내놓으면 director를 그만큼 돌려야 한다. 회전행렬은 로드리게스(Rodrigues) 공식으로 만든다.

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

$\Theta$ 는 $\boldsymbol{\theta}$ 의 반대칭 행렬, $\theta = |\boldsymbol{\theta}|$ 다. 증분이 작다고 $R \approx I + \Theta$ 로 끊는 코드가 흔하다. 이 행렬은 직교행렬이 아니다. $\det(I + \Theta) = 1 + \theta^{2} \neq 1$ 이므로 매 스텝 director가 조금씩 길어진다.

아래 시뮬레이션에서 직접 조작해보자.

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.

d.theta를 0.10에 두고 진행시키면 빨간 1차 절단 화살표가 점선 단위원 밖으로 나선을 그리며 빠져나간다. 초록 닫힌 형태는 원 위에 정확히 남는다. d.theta를 절반으로 줄이면 드리프트도 절반이 된다. 즉 이 오차는 증분 크기에 1차다.

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 는 급수 절단, 0 은 닫힌 형태"""
    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):
    """중립면 법선을 y축 둘레로 스텝당 dth 라디안씩 돌린다"""
    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

1차 절단은 30스텝 만에 director를 16.1% 늘였다. 항 하나를 더 붙인 2차는 0.04%로 떨어진다. 길이가 400배 정확해졌는데 각도 오차는 오히려 2차 쪽이 크다. 두 오차는 별개다.

4번 표가 더 중요하다. 증분을 10분의 1로 줄여도 오차는 10분의 1로만 준다. 하중 스텝을 잘게 쪼개는 것으로는 이 문제를 못 없앤다. 닫힌 형태를 쓰거나, 매 스텝 director를 정규화하거나, 회전을 쿼터니언으로 들고 다녀야 한다.

$S$ 인지 $\sigma$ 인지 코드에서 가려내는 법#

남의 대변형 코드를 열었을 때 응력 변수의 정체는 이름으로는 알 수 없다. 세 곳을 보면 갈린다.

응력이 강성행렬로 들어가는 자리를 본다. $B$ 행렬이 $\partial E / \partial u$ 로 만들어져 있으면 곱해지는 응력은 $S$ 다. $B$ 가 $\partial \varepsilon / \partial u$ 면 $\sigma$ 다. 구성행렬 변환에서 봤듯 $B$ 의 정의가 응력의 정체를 결정한다.

적분 시 야코비안을 본다. $\int_{V_0} \cdots \, dV_0$ 처럼 기준 배치 부피로 적분하면서 $J$ 를 곱하지 않으면 $S$ 다. $J$ 를 곱하고 있다면 $\sigma$ 를 기준 배치로 되돌리는 중이다.

출력 루틴을 본다. 후처리에서 von Mises를 계산하기 직전에 $\sigma = J^{-1}FSF^{T}$ 변환이 들어 있으면, 내부적으로는 $S$ 로 굴러가고 있었다는 뜻이다. 이 변환 없이 $S$ 의 성분을 그대로 von Mises에 넣어 그리는 코드가 있다. 작은 변형에서는 티가 안 나고 늘임이 20%를 넘으면 갈라진다.

셋 다 확인이 안 되면 테스트가 남는다. 요소 하나를 강체 회전만 시키고 응력이 0으로 유지되는지 본다. 30°에서 13%가 나오면 어디선가 $F$ 를 제곱하지 않았다는 뜻이다.