요소를 45° 비틀었더니 강성이 2.5배가 됐다 — 메트릭 텐서와 구성 텐서 좌표변환

직교 정규 기저에서는 공변 성분과 반변 성분이 같은 숫자다. 그 우연이 깨지는 순간, 구성행렬을 그대로 쓴 코드가 틀리기 시작한다.

같은 변형을 두 번 재는데 답이 다르게 나오는 상황이 있다. 요소는 그대로, 재료도 그대로, 실제로 일어난 변형도 그대로다. 바꾼 것은 그 변형을 적어 두는 좌표계뿐이다. 그런데 변형에너지가 2.5배로 늘어난다. 이 글은 그 2.5배가 어디에서 오는지, 공변·반변 기저와 메트릭 텐서로 짚고, 4차 구성 텐서를 어떻게 변환해야 숫자가 제자리로 돌아오는지 Python으로 확인한다.

변형은 그대로인데 에너지가 2.5배가 됐다#

쉘 요소의 강성행렬은 대개 이렇게 조립된다. 가우스 적분점에서 변형률-변위 행렬 $\mathbf{B}$ 를 만들고, 구성행렬 $\mathbf{C}$ 를 곱하고, $\mathbf{B}^{\mathsf T}\mathbf{C}\mathbf{B}$ 를 적분한다. 문제는 이 두 행렬이 서로 다른 좌표계에서 태어난다는 데 있다.

$\mathbf{C}$ 는 재료의 성질이다. 그래서 물성 시험이 이루어진 좌표계, 즉 쉘 표면에 붙인 국소 직교 데카르트 좌표계에서 정의된다. 반면 $\mathbf{B}$ 는 형상함수의 미분에서 나온다. 형상함수는 요소의 자연좌표계 $(r^1, r^2, r^3)$ 에서 쓰여 있고, 요소가 휘거나 찌그러지면 이 좌표계는 직교하지도, 정규화되지도 않는다.

두 좌표계가 우연히 일치하는 경우가 있다. 요소가 직사각형이고 표면이 평평할 때다. 패치 테스트를 평면 모델로만 돌리면 이 우연 위에서만 검증하게 된다. 곡면에 얹는 순간 우연은 깨진다.

기저가 직교하지 않으면 성분이 두 벌 생긴다#

자연좌표 $r^i$ 가 변할 때 물리 공간의 위치 $\mathbf{x}$ 가 움직이는 방향, 그것이 공변 기저(covariant basis)다.

\mathbf{g}_i = \frac{\partial \mathbf{x}}{\partial r^i}

$\mathbf{x}$ 는 데카르트 위치벡터, $r^i$ 는 요소의 자연좌표다. 이 세 벡터는 서로 직교하지 않고 길이도 1이 아니다.

직교하지 않는 기저 위에서 벡터 하나를 성분으로 적는 방법은 둘이다. 첫째, 기저를 따라 평행사변형으로 분해한다. 이때 나오는 계수가 반변 성분(contravariant component) $v^i$ 다. 둘째, 각 기저 벡터 위로 수직으로 떨어뜨린다. 그 사영이 공변 성분(covariant component) $v_i$ 다.

\mathbf{v} = v^i \mathbf{g}_i, \qquad v_i = \mathbf{v}\cdot\mathbf{g}_i

직교 정규 기저에서는 두 구성이 같은 점에 도착한다. 그래서 데카르트 좌표계만 쓰던 사람에게는 이 구분이 아예 보이지 않는다. 아래 시뮬레이션에서 직접 조작해보자.

skew of g_2 [deg]30| g_2 |1.40

v angle [deg]35

The white arrow is one physical vector. Blue dashes slide it along the green base vectors (contravariant vⁱ); yellow dashes drop it perpendicular onto them (covariant v_i). Push the skew slider to 0 and |g₂| to 1.00: the two rows of numbers merge and max | v^i - v_i | turns green. Anywhere else they differ, and only the metric row keeps the length right.

skew를 0, |g_2|를 1.00에 두면 파란 평행사변형과 노란 수선이 같은 점에서 만나고 max | v^i - v_i |가 초록으로 바뀐다. 슬라이더를 조금만 움직여도 두 줄의 숫자가 갈라진다. 흰 화살표는 내내 같은 벡터라는 점을 보면서 조작하는 것이 핵심이다.

메트릭 텐서가 잃어버린 길이를 돌려준다#

두 성분을 잇는 것이 메트릭 텐서다.

g_{ij} = \mathbf{g}_i \cdot \mathbf{g}_j, \qquad v_i = g_{ij}\, v^j

$g_{ij}$ 는 기저 벡터들의 길이와 사잇각을 한 행렬에 담은 것이다. 대각 성분 $g_{11}, g_{22}, g_{33}$ 은 각 기저의 길이 제곱이고, 비대각 성분은 사잇각의 코사인에 길이를 곱한 값이다.

길이를 재려면 반드시 이 행렬을 거쳐야 한다.

|\mathbf{v}|^2 = g_{ij}\, v^i v^j = v_i v^i

두 번째 등식이 중요하다. 공변 성분과 반변 성분을 짝지어 곱하면 메트릭이 저절로 사라진다. 반면 반변 성분끼리 제곱해서 더하면( $v^i v^i$ ) 그것은 길이가 아니다. viz의 맨 아랫줄이 그 값을 빨갛게 보여준다.

연속체 역학이 응력에 반변, 변형률에 공변을 쓰는 이유가 여기 있다. 가상일은 스칼라라야 하고, 2차 Piola-Kirchhoff 응력의 반변 성분 $S^{ij}$ 와 Green-Lagrange 변형률의 공변 성분 $E_{ij}$ 를 짝지어야 $S^{ij}E_{ij}$ 가 좌표계와 무관해진다. 격자 메트릭이 유동 해석에서 하는 일과 같다 — 곡선 좌표 변환과 격자 메트릭에서 다룬 야코비안이 여기서는 기저 벡터의 정체다.

반변 기저는 야코비안 역행렬의 행이다#

공변 성분을 반변으로 되돌리려면 $g^{ij} = [g_{ij}]^{-1}$ 이 필요하고, 이 행렬로 만든 기저가 반변 기저다.

\mathbf{g}^i = g^{ij}\mathbf{g}_j, \qquad \mathbf{g}^i \cdot \mathbf{g}_j = \delta^i_j

$\delta^i_j$ 는 크로네커 델타다. 즉 $\mathbf{g}^1$ 은 $\mathbf{g}_2, \mathbf{g}_3$ 둘 다에 수직이면서 $\mathbf{g}_1$ 과의 내적이 1이 되도록 길이가 맞춰진 벡터다.

구현에서는 역행렬 두 번을 돌릴 필요가 없다. 공변 기저를 열로 쌓은 행렬을 $\mathbf{J}$ 라 하면,

\mathbf{J} = [\,\mathbf{g}_1\;\mathbf{g}_2\;\mathbf{g}_3\,], \qquad \mathbf{g}^i = \big(\mathbf{J}^{-1}\big)_{i\bullet}

$\big(\mathbf{J}^{-1}\big)_{i\bullet}$ 는 $\mathbf{J}^{-1}$ 의 $i$ 번째 행이다.

$\mathbf{J}$ 는 자연좌표에서 물리좌표로 가는 야코비안 그 자체다. 이미 적분점마다 $\det\mathbf{J}$ 를 구하려고 계산해 둔 행렬이다. 반변 기저는 덤으로 나온다.

4차 텐서는 방향코사인을 네 번 곱한다#

이제 본론이다. 국소 데카르트 기저 $\mathbf{e}_p$ 에서 정의된 구성 텐서 $C^{pqrs}$ 를 자연좌표계로 옮긴다. 2차 텐서를 옮길 때 방향코사인을 두 번 곱했다면, 4차 텐서는 네 번 곱한다.

\tilde{C}^{ijkl} = (\mathbf{g}^i\!\cdot\mathbf{e}_p)(\mathbf{g}^j\!\cdot\mathbf{e}_q)(\mathbf{g}^k\!\cdot\mathbf{e}_r)(\mathbf{g}^l\!\cdot\mathbf{e}_s)\, C^{pqrs}

인덱스 $i,j,k,l$ 은 자연좌표, $p,q,r,s$ 는 국소 데카르트 좌표를 가리킨다. 변형률 쪽은 방향이 반대다.

\tilde{\varepsilon}_{ij} = (\mathbf{g}_i\!\cdot\mathbf{e}_p)(\mathbf{g}_j\!\cdot\mathbf{e}_q)\,\varepsilon_{pq}

구성 텐서에는 $\mathbf{J}^{-1}$ 이, 변형률에는 $\mathbf{J}$ 가 붙는다. 그래서 둘을 축약하면 $\mathbf{J}$ 들이 정확히 상쇄되고 에너지가 불변으로 남는다. 반대로 말하면, 변형률만 변환하고 구성 텐서를 그대로 두면 $\mathbf{J}$ 가 네 개 남는다. 그 네 개가 아래에서 볼 오차의 정체다.

skew of g_2 [deg]30| g_2 |1.00

The blue element on the left never changes its physics — only the coordinate lines you measure it with do. Drag the skew slider away from 0, or push |g₂| off 1.00, and watch the red bar climb while the green one stays put. Bring both back to the orthonormal setting and the two bars lock together: that is the only configuration where skipping the transform is free.

skew와 |g_2|를 움직이면 왼쪽 요소의 변형 모양은 그대로인데 오른쪽 빨간 막대만 자란다. stretch/shear/mixed를 바꿔가며 위쪽 오차 곡선의 모양이 어떻게 달라지는지 보는 것이 관찰 포인트다. 전단 모드에서 오차가 가장 빨리 커진다.

Python으로 재본 뒤틀림 각도별 에너지#

평면응력 등방 재료에 변형률 하나를 고정해 놓고, 좌표계만 비틀면서 변형에너지를 두 방식으로 계산했다. 원본 문서와 같은 물성을 썼다( $E = 2.1\times10^6$ , $\nu = 0.3$ ).

import numpy as np
 
def natural_basis(skew_deg, stretch=1.0):
    """열(column)이 공변 기저 g_i = dx/dr^i 인 야코비안 J"""
    a = np.deg2rad(skew_deg)
    g1 = np.array([1.0, 0.0])
    g2 = stretch * np.array([np.sin(a), np.cos(a)])
    return np.column_stack([g1, g2])
 
def plane_stress_tensor(E=2.1e6, nu=0.3):
    """국소 데카르트 좌표계의 등방 평면응력 4차 텐서 C^{pqrs}"""
    lam = E * nu / (1.0 - nu**2)
    mu = E / (2.0 * (1.0 + nu))
    d = np.eye(2)
    return (lam * np.einsum('pq,rs->pqrs', d, d)
            + mu * (np.einsum('pr,qs->pqrs', d, d) + np.einsum('ps,qr->pqrs', d, d)))
 
def rotate_fourth_order(C, Jinv):
    """C^{ijkl} = (g^i.e_p)(g^j.e_q)(g^k.e_r)(g^l.e_s) C^{pqrs}, g^i = J^-1 의 i번째 행"""
    return np.einsum('ip,jq,kr,ls,pqrs->ijkl', Jinv, Jinv, Jinv, Jinv, C)
 
def strain_energy(C, eps):
    return 0.5 * np.einsum('pqrs,pq,rs->', C, eps, eps)
 
def to_voigt2d(C):
    """2D Voigt: (00,11,01) -> 3x3"""
    idx = [(0, 0), (1, 1), (0, 1)]
    return np.array([[C[p, q, r, s] for (r, s) in idx] for (p, q) in idx])
 
# 물리적 변형률 상태 하나 (국소 데카르트 성분). 좌표계를 어떻게 잡든 이 텐서는 그대로다.
eps_cart = np.array([[1.0e-3, 4.0e-4],
                     [4.0e-4, -6.0e-4]])
C_cart = plane_stress_tensor()
U_ref = strain_energy(C_cart, eps_cart)
 
print("skew   g11     g12     g22    | U_correct     U_naive      err%")
print("-" * 68)
for skew in [0, 5, 10, 15, 20, 30, 40, 45]:
    J = natural_basis(skew)
    g = J.T @ J                       # 메트릭 텐서 g_ij
    Jinv = np.linalg.inv(J)           # 행 = 반변 기저 g^i
    eps_nat = J.T @ eps_cart @ J      # 공변 변형률 성분
    C_nat = rotate_fourth_order(C_cart, Jinv)
    U_ok = strain_energy(C_nat, eps_nat)
    U_bad = strain_energy(C_cart, eps_nat)   # 변환을 빠뜨린 코드
    err = 100.0 * (U_bad - U_ok) / U_ok
    print(f"{skew:3d}   {g[0,0]:.3f}  {g[0,1]:+.3f}  {g[1,1]:.3f} |"
          f" {U_ok:.6e}  {U_bad:.6e}  {err:+8.2f}")
 
print()
print("orthogonal (skew=0), only |g2| stretched")
for st in [1.0, 1.5, 2.0]:
    J = natural_basis(0, stretch=st)
    Jinv = np.linalg.inv(J)
    eps_nat = J.T @ eps_cart @ J
    U_ok = strain_energy(rotate_fourth_order(C_cart, Jinv), eps_nat)
    U_bad = strain_energy(C_cart, eps_nat)
    print(f"  stretch={st:.1f}  g22={(J.T@J)[1,1]:.2f}  err% = {100*(U_bad-U_ok)/U_ok:+9.2f}")
 
print()
print(f"Cartesian reference U_ref   = {U_ref:.6e}")
J = natural_basis(30)
Jinv = np.linalg.inv(J)
C_nat = rotate_fourth_order(C_cart, Jinv)
eps_nat = J.T @ eps_cart @ J
print(f"skew=30 after transform     = {strain_energy(C_nat, eps_nat):.6e}  (invariant)")
 
# Voigt로 접어도 같은 값이 나오는지
Cv = to_voigt2d(C_nat)
ev = np.array([eps_nat[0, 0], eps_nat[1, 1], 2.0 * eps_nat[0, 1]])
print(f"skew=30 via Voigt 3x3       = {0.5 * ev @ Cv @ ev:.6e}")
 
# Voigt 공간에서 변형률 변환행렬 A: e_v(nat) = A e_v(cart)
def voigt_map(J):
    cols = []
    for e in (np.array([[1, 0], [0, 0]]), np.array([[0, 0], [0, 1]]), np.array([[0, .5], [.5, 0]])):
        n = J.T @ e @ J
        cols.append([n[0, 0], n[1, 1], 2 * n[0, 1]])
    return np.array(cols).T
 
A = voigt_map(J)
Cv_cart = to_voigt2d(C_cart)
Ai = np.linalg.inv(A)
print("Voigt congruence C_nat = A^-T C_cart A^-1  residual =",
      f"{np.max(np.abs(Ai.T @ Cv_cart @ Ai - Cv)):.3e}")

skew   g11     g12     g22    | U_correct     U_naive      err%
--------------------------------------------------------------------
  0   1.000  +0.000  1.000 | 1.412308e+00  1.412308e+00     +0.00
  5   1.000  +0.087  1.000 | 1.412308e+00  1.486003e+00     +5.22
 10   1.000  +0.174  1.000 | 1.412308e+00  1.585621e+00    +12.27
 15   1.000  +0.259  1.000 | 1.412308e+00  1.722495e+00    +21.96
 20   1.000  +0.342  1.000 | 1.412308e+00  1.906550e+00    +35.00
 30   1.000  +0.500  1.000 | 1.412308e+00  2.437219e+00    +72.57
 40   1.000  +0.643  1.000 | 1.412308e+00  3.163176e+00   +123.97
 45   1.000  +0.707  1.000 | 1.412308e+00  3.567692e+00   +152.61
 
orthogonal (skew=0), only |g2| stretched
  stretch=1.0  g22=1.00  err% =     +0.00
  stretch=1.5  g22=2.25  err% =   +105.60
  stretch=2.0  g22=4.00  err% =   +407.84
 
Cartesian reference U_ref   = 1.412308e+00
skew=30 after transform     = 1.412308e+00  (invariant)
skew=30 via Voigt 3x3       = 1.412308e+00
Voigt congruence C_nat = A^-T C_cart A^-1  residual = 9.313e-10

읽을 것이 셋 있다.

첫째, skew=0 행에서 오차가 정확히 0이다. 직사각형 요소만으로 돌린 검증은 이 버그를 절대 잡지 못한다. 둘째, 15° 뒤틀림에서 벌써 22%다. 실제 곡면 격자에서 흔한 각도다. 셋째, 직교인데도 $|\mathbf{g}_2|$ 만 1.5배로 늘리면 오차가 105%다. 뒤틀림이 아니라 메트릭이 단위행렬이 아닌 것이 원인이라는 뜻이다.

Voigt로 접으면 6×6 행렬 하나가 된다#

인덱스 네 개짜리 배열을 코드에 그대로 들고 다니지는 않는다. 응력·변형률 텐서가 대칭이므로 독립 성분은 6개뿐이고, 4차 텐서는 $6\times6$ 행렬로 접힌다(위 코드는 2D라 $3\times3$ ).

접은 뒤에도 변환은 살아 있다. 변형률 Voigt 벡터가 $\tilde{\boldsymbol\varepsilon}_v = \mathbf{A}\,\boldsymbol\varepsilon_v$ 로 변한다면, 에너지가 불변이어야 하므로 구성행렬은 합동변환을 받는다.

\tilde{\mathbf{C}}_v = \mathbf{A}^{-\mathsf T}\,\mathbf{C}_v\,\mathbf{A}^{-1}

$\mathbf{A}$ 의 성분은 방향코사인의 곱이다. 위 출력 마지막 줄의 잔차 $9.3\times10^{-10}$ 이 이 식이 4차 텐서 축약과 같은 답을 준다는 확인이다. 실제 쉘 코드에서 흔히 보는 T 행렬이 바로 이것이고, 전단 보정 계수 5/6과 평면응력 가정( $\sigma_{33}=0$ )은 $\mathbf{C}_v$ 를 만들 때 국소 데카르트 좌표계에서 먼저 반영한다. 순서를 바꾸면 안 된다 — 평면응력 조건은 두께 방향이 정의된 그 좌표계에서만 의미가 있다.

유한체적법에서 같은 실수가 나오는 자리#

이 실수는 구조해석 코드만의 것이 아니다. 곡선격자 유한체적법에서 점성응력 텐서를 계산할 때 같은 구조가 나온다. 변형률률(strain rate) 텐서를 자연좌표 미분으로 얻어 놓고, 뉴턴 점성 법칙 $\tau = 2\mu \mathbf{D}$ 를 데카르트 형태 그대로 적용하면 정확히 위 표의 오차가 생긴다.

셋을 확인하면 된다. 첫째, 지금 들고 있는 텐서 성분이 물리 성분인가 공변/반변 성분인가. 둘째, 축약할 때 위 인덱스와 아래 인덱스가 짝을 이루는가. 셋째, 검증 케이스에 뒤틀린 요소가 하나라도 있는가.

세 번째가 실질적으로 가장 중요하다. 비직교 확산 플럭스 보정에서도 그랬듯, 비직교성이 만드는 오차는 직교 격자 검증을 100% 통과한 뒤에 나타난다. 쉘 요소의 전단 잠김과 MITC 타잉이 요소의 공식을 고치는 이야기였다면, 이 글은 그 공식을 어느 좌표계에서 읽느냐의 이야기다. 둘 다 틀려도 평면 패치 테스트는 통과한다.