FESADev/docs/shell-transient-heat-transfer/formulation.md

# Shell 과도 열전달 정식화

- 상태: **2026-09-29 승인된 물리·방법을 반영한 구현용 정식화 / 독립 numerical review 대기**.
- 담당 범위: `FEM-REQ-SHELL-THERMAL-001`–`018`의 수식·물리 의미와 검증 함의.
- 기준: [Feature PRD](PRD.md), [Feature ADR](ADR.md), [프로젝트 지침](../../AGENTS.md),
  [프로젝트 PRD](../PRD.md), [공통 ADR](../ADR.md#adr-028), [ARCHITECTURE](../ARCHITECTURE.md).
- 아래 수식은 승인된 FESA 모델이다. Abaqus 내부 알고리즘과 동등하다는 주장이 아니다.
  이전 Proposed 보고서의 검산은 제한된 과거 evidence이며 현재 solver 검증을 대신하지 않는다.

## 1. 적용 가정과 차원

고정된 3차원 midsurface의 4절점 Q4 요소에 두께 방향 온도점 5개를 둔다.
요소마다 유한한 두께 `h > 0`, 등방성 열전도율 `k > 0`, 밀도 `ρ > 0`, 비열 `c > 0`가
상수이며 centered 단면 전체는 같은 재료다. 요소 집합별 재료·두께를 배정할 수 있다.
단일 uncoupled transient step만 지원한다. 변형·연성·대류·복사·내부 발열·적층·분포
표면 열유속·분포 edge 입력·시간 변화 입력·다중 step은 제외한다.

일관된 입력 단위계를 사용한다. SI 예시는 `h [m]`, `T [K]`, `t [s]`,
`k [W/(m·K)]`, `ρ [kg/m³]`, `c [J/(kg·K)]`, 열유속 `q [W/m²]`이다.
상수 물성에서는 온도 차에 같은 척도를 쓰는 섭씨도 가능하다. 숫자로 단위계를 추정하지 않는다.
온도 기준 이동에 불변인 에너지 증분을 사용하며 절대 내부에너지의 영점을 추정하지 않는다.

기하학적으로 얇은 shell의 product-domain 근사로 면내·두께 전도를 모두 유지한다.
여기서 ‘얇음’은 **두께 방향 등온 가정이 아니다**. 정적인 geometry이므로 구조 회전 DOF,
기하 비선형 tangent, 재료 history update는 없다. 유일한 진화 state는 온도다.

## 2. Geometry, 방향과 접속의 의미

자연좌표 `(ξ,η) ∈ [-1,1]²`, 절점 순서 `(-1,-1),(1,-1),(1,1),(-1,1)`를 사용한다.
절점 좌표를 `x_i`, midsurface를 `r=Σ_i N_i x_i`로 두고 다음을 계산한다.

\[
a_1=r_{,\xi},\quad a_2=r_{,\eta},\quad
g_{\alpha\beta}=a_\alpha\cdot a_\beta,\quad
J_s=\|a_1\times a_2\|=\sqrt{\det g},\quad
n=(a_1\times a_2)/J_s,\quad
a^\alpha=\sum_\beta(g^{-1})_{\alpha\beta}a_\beta.
\]

`α,β=1,2`; `a^α`는 reciprocal tangent이다. 실제 평가점마다 metric을 계산하며 warped
Q4를 요소 중심의 평면에 투영하지 않는다. 출력 직교축은 `e1=a1/|a1|`, `e2=n×e1`, `e3=n`이다.
두께좌표 `ζ ∈ [-1,1]`, `z=hζ/2`; 음의 normal 쪽이 bottom, 양의 쪽이 top이다.
표시 위치는 `x=r+zn`으로 정의하되, 연산자는 아래 근사를 쓴다.

\[
\nabla_* T=a^1T_{,\xi}+a^2T_{,\eta}+n\frac{2}{h}T_{,\zeta},\qquad
dV_* = J_s\frac h2\,d\xi\,d\eta\,d\zeta.
\]

정확한 offset-volume의 접선 `aα+z n_,α`, 두께별 metric·면적 변화는 생략한다.
따라서 `h |κ1|`, `h |κ2|`가 작고 열장이 변하는 면내 길이에 비해 두께가 작은 영역을
대상으로 하는 물리 근사다. 이는 승인되지 않은 curvature/thickness cutoff를 입력에 추가하는
조건이 아니다. 평판·warped patch·곡면 mesh·fold를 지원하며 적분 차수를 높여도 생략한
offset-volume 모델링 오차는 없어지지 않는다.

Geometry의 필수 조건은 유한 좌표, 서로 다른 4절점, `J_s>0`, 비특이 metric, 자기 교차 없는
유효 patch와 일관된 orientable connectivity이다. 적분점 검사만으로 전체 patch 유효성을
증명하지 않는다. 아래 whole-patch certificate를 사용하며 작은 양수 Jacobian을 clamp하지 않는다.

**승인된 접속:** 인접 요소는 같은 source node의 같은 `ζ` 온도 DOF를 공유한다.
매끈한 면에서 이는 온도 연속성과 약형 열률 평형을 준다. 접힌 판에서는 요소별 normal을
유지하고 공유 edge의 `ζ`별 온도만 같게 묶는 이상적 무저항·무용량 접합이다.
서로 다른 두께에서는 같은 normalized `ζ`를 연결하므로 실제 `z` 좌표가 같다는 뜻이 아니다.
실제 step 단차의 노출면, 유한 접합부 부피, 국소 spreading resistance는 해상하지 않는다.
이는 승인된 이상화이며 상세 3D 접합부 온도장·열유속의 재현을 의미하지 않는다.
공유 edge의 요소 순서는 서로 반대여야 하고 bottom/top의 연속된 의미를 유지해야 한다.
온도는 scalar이므로 normal 사이의 회전행렬을 곱하지 않는다. 법선 반전은 해당 두께 station
`a ↔ 6-a`와 면 하중/경계조건을 함께 바꿔야 같은 물리 문제가 된다.

### 2.1 Whole-patch geometry certificate

먼저 네 source node identity와 좌표가 각각 서로 다르고 유한함을 확인한다. Bilinear map을
`r=a+bξ+cη+dξη`로 쓰면 다음은 근사가 아닌 항등식이다.

\[
b=(-x_1+x_2+x_3-x_4)/4,\quad
c=(-x_1-x_2+x_3+x_4)/4,\quad
d=(x_1-x_2+x_3-x_4)/4,
\]
\[
w(\xi,\eta)=a_1\times a_2=w_0+u\xi+v\eta,\qquad
w_0=b\times c,\quad u=b\times d,\quad v=d\times c.
\]

따라서 `J_s²=||w||²`의 사각형 최소값은 convex quadratic minimization이다.
정확산술에서 다음 유한 후보 집합으로 최소점을 구할 수 있다.

1. 네 모서리 `(±1,±1)`.
2. `ξ=s=±1`의 edge: `v·v>0`이면
   `η=clip(-v·(w0+su)/(v·v),[-1,1])`; `v=0`이면 모서리로 충분하다.
   `η=s`의 edge는 `u,v`와 `ξ,η`를 교환한다. 이 clip은 constrained minimizer의
   수학적 projection이며 geometry/Jacobian의 clamp가 아니다.
3. `D=(u·u)(v·v)-(u·v)²>0`이면
   `[[u·u,u·v],[u·v,v·v]] [ξ,η]ᵀ=-[u·w0,v·w0]ᵀ`의 해가 사각형 안에 있을 때 추가한다.
   `D=0`이면 minimum의 affine level set이 경계에 도달하므로 edge 후보로 충분하다.

최소점 `p*`와 `w*=w(p*)`에 대해 convex projection 성질은
`w*·w(p)≥||w*||²`이다. 그러므로 **어떤 비영 벡터 `s`에 대해서든 네 모서리 모두
`s·w_corner>0`를 증명하면 전체 patch의 `J_s>0`가 증명된다.** 역으로 유효 patch는
정확산술에서 `s=w*`가 이 certificate를 준다. 고정 normal에 대한 양의 projection을
가정하는 방법이 아니라 각 patch에서 certificate 방향을 찾는 방법이다.

자기 교차도 함께 배제된다. 임의의 두 자연좌표 `p,q`에 대해 bilinear 항등식
`r(p)-r(q)=Dr((p+q)/2)(p-q)`가 성립한다. 중점은 사각형 안에 있고 `J_s>0`이면
`Dr`의 rank는 2이므로 `r(p)=r(q)`는 `p=q`를 뜻한다. `g`는 Gram 행렬이므로
동시에 positive definite다. 별도 corner-only 또는 quadrature-only 판정은 충분하지 않다.

구현의 부동소수점 판정은 다음처럼 certificate와 후보 탐색을 분리한다.

- IEEE binary64 round-to-nearest, gradual underflow를 전제로 한다. Fast-math,
  flush-to-zero 또는 검증하지 않은 연산 재결합을 사용하지 않는다. 좌표 차·곱·합을
  powers-of-two로 재척도하여 overflow/underflow를 피한다. 길이 척도는 실제 좌표
  차이에서 얻으며 `max(1,scale)` 같은 단위 의존 floor를 두지 않는다.
- 각 입력 double을 singleton interval로 보고 `+,-,×,÷` 각각의 계산 결과를
  `nextafter(result,-∞)`와 `nextafter(result,+∞)`로 바깥으로 넓힌다.
  곱/나눗셈은 endpoint 조합의 extrema를 사용하며 0을 포함하는 divisor는 사용하지
  않는다. Exact powers-of-two 변환도 subnormal 손실이 있으면 enclosure에 포함한다.
  이 규칙은 cancellation을 포함해 연산 횟수에 맞는 bound를 전파한다.
- 위 후보를 통상 double로 구해 `s=w(candidate)` 방향을 시도한다. 후보의 작은
  오차나 사각형 바깥 위치는 certificate의 건전성을 해치지 않는다. 후보 탐색의 `D`
  또는 division이 불확실하면 유한한 다른 후보를 시도한다. 각 `s` 자체는 고정 double
  벡터로 취급하고, 원래 좌표로부터 interval로 다시 구한 네 `w_corner`에 대해
  `lower(s·w_corner)>0`인지 검증한다. 최종 판정은 후보 최소값의 부호에 의존하지 않는다.
- 모두 양수이면 `J_s≥min_corner lower(s·w_corner)/upper(||s||)>0`라는
  whole-patch lower bound가 생긴다. Certificate가 없으면 invalid/uncertain geometry로
  fail closed하고 source element와 실패 bound를 진단한다. 임의 aspect/warp cutoff,
  고정 epsilon margin 또는 arbitrary-precision dependency는 필요하지 않다.
  수치적으로 증명 불가능한 경우를 정확한 기하학적 퇴화라고 단정하지 않는다.

평가점에서도 유한한 metric·gradient를 확인한다. Certificate 통과가 극단적 물성·두께나
metric 역행렬의 유한한 계산을 보장하지는 않으며, 그 실패를 normal 수정으로 숨기지 않는다.

### 2.2 Mesh topology와 source 접속

Source edge key는 정렬한 두 node identity이며 각 incidence에는 원래 directed pair를
보존한다. Incidence는 1 또는 2이고, 2이면 방향이 반대여야 한다. 3개 이상의 요소가
공유하는 edge, 같은 방향의 공유 edge, 중복 element patch는 거부한다. 자동 normal flip이나
좌표 기반 node welding은 하지 않는다. Disconnected components는 각각 유효하다.

각 vertex의 link graph를 검사한다. Link vertex는 해당 node에 닿는 source edge이고,
각 incident Q4는 그 두 edge 사이에 link edge 하나를 만든다. Link는 내부 vertex에서
하나의 cycle, boundary vertex에서 하나의 path여야 한다. 끊어진 link, 분기, 중복 연결은
vertex-only fan/nonmanifold junction이므로 거부한다. Fold의 각도 자체는 거부 이유가 아니다.
T-junction은 source incidence만으로 찾을 수 없으므로 mesh node가 다른 edge의 열린
선분 내부에 있는지 검사한다: `x=p+λ(q-p), 0<λ<1`의 collinearity 및 interval을
scale-aware bound로 판정하고 걸쳐 있는 경우 uncertain 접속으로 진단한다. 접속을 새로
생성하지 않는다. 같은 edge를 부분적으로 덮는 중복 선분도 이 endpoint 검사에 포함한다.
이 검사는 patch 내부 injectivity와 별개이며 mesh 사이 접촉을 열접속으로 추정하지 않는다.

## 3. 20 DOF와 온도 보간

절점별 온도점은 `ζ_a=(-1,-1/2,0,1/2,1)`, `a=1,…,5`이다.
로컬 순서는 `p=5(i-1)+a`, 즉 node-major 20 DOF다. 외부 DOF는 `10+a`로 대응한다.
따라서 `11=bottom`, `13=midsurface`, `15=top`이다. 기하 절점 위치와 DOF11의 물리 위치를
혼동하지 않는다. 외부 위치 의미는 승인된 DS4 입력 계약을 따른다.

면내 shape와 미분은 각 절점의 부호 `(ξ_i,η_i)`를 사용한다.

\[
N_i=\tfrac14(1+\xi_i\xi)(1+\eta_i\eta),\quad
N_{i,\xi}=\tfrac14\xi_i(1+\eta_i\eta),\quad
N_{i,\eta}=\tfrac14\eta_i(1+\xi_i\xi).
\]

두께는 `I−=[-1,0]`, `I+=[0,1]`의 **두 C0 quadratic 구간**을 사용한다.
단일 quartic보다 각 구간의 저차 보간·국소 미분을 명료하게 유지하며, 중간 온도 1개를 공유한다.
구간 좌표는 `s=2ζ+1` 또는 `s=2ζ−1`; 두 경우 `s∈[-1,1]`, `ds/dζ=2`이다.

\[
\ell_1(s)=\tfrac12s(s-1),\quad \ell_2(s)=1-s^2,\quad
\ell_3(s)=\tfrac12s(s+1),\qquad
\ell'_1=s-\tfrac12,\quad\ell'_2=-2s,\quad\ell'_3=s+\tfrac12.
\]

`I−`에서는 `(L1,L2,L3)=(ℓ1,ℓ2,ℓ3)`, `L4=L5=0`;
`I+`에서는 `(L3,L4,L5)=(ℓ1,ℓ2,ℓ3)`, `L1=L2=0`이다.
활성 함수의 `L_,ζ=2ℓ'`; `ζ=0`의 함수값은 유일하지만 미분은 좌·우로 구분한다.

\[
H_{ia}=N_iL_a,\qquad T_h=\sum_{i,a}H_{ia}T_{ia}=\mathbf H\mathbf T_e,\qquad
B_{ia}=L_a(a^1N_{i,\xi}+a^2N_{i,\eta})+n\frac2hN_iL_{a,\zeta}.
\]

`H`는 `1×20`, `B`는 global Cartesian `3×20`, `∇_*T_h=B T_e`다.
`ΣH=1`, `ΣB=0`; Kronecker interpolation과 내부 station 온도 연속성을 확인한다.
두 구간은 수치 보간 구간이며 물성 layer나 물리 계면이 아니다.

## 4. 지배식, 약형과 20×20 행렬

Fourier 법칙 `q=−k∇_*T`, 저장식 `ρc Tdot + div_* q=0`를 사용한다.
외향 열유속 `q·ν`와 반대인 `q_in=−q·ν`를 양의 입력으로 정의한다.
자유 test function `v`에 대해 product-domain 약형은 다음과 같다.

\[
\int_{V_*}\rho c\,v\dot T\,dV_*+
\int_{V_*}k\nabla_*v\cdot\nabla_*T\,dV_*
=\sum_j v_jP_j.
\]

지정 온도 외의 하중 없는 경계는 자연 단열이다. 집중 입력 `P_j`는 승인된 bottom/top
절점의 이산 heat rate이며 연속체 point-source 근방의 유한 heat flux를 보장하지 않는다.

\[
\mathbf C_e=\int\rho c\,\mathbf H^T\mathbf H\,J_s\frac h2\,d\xi d\eta d\zeta,
\qquad
\mathbf K_e=\int k\,\mathbf B^T\mathbf B\,J_s\frac h2\,d\xi d\eta d\zeta.
\]

`C_e [J/K]`는 consistent 열용량, `K_e [W/K]`는 열전도 행렬이다.
`C_e` 대칭 양정치, `K_e` 대칭 양반정치와 유일한 요소 상수온도 null mode를 기대한다.
Reduced integration, hourglass stiffness, capacity lumping, 인공 확산은 사용하지 않는다.

명시적인 component식과 분리 계산에 필요한 작은 행렬은 다음과 같다.

\[
M_{ij}=\int_{[-1,1]^2}N_iN_jJ_s\,d\xi d\eta,\quad
S_{ij}=\int_{[-1,1]^2}\sum_{\alpha\beta}N_{i,\alpha}g^{\alpha\beta}N_{j,\beta}J_s\,d\xi d\eta,
\quad Z^0_{ab}=\int_{-1}^1L_aL_b\,d\zeta,\quad Z^1_{ab}=\int_{-1}^1L_{a,\zeta}L_{b,\zeta}\,d\zeta,
\]
\[
(C_e)_{ia,jb}=\rho c\frac h2M_{ij}Z^0_{ab},\qquad
(K_e)_{ia,jb}=k\left(\frac h2S_{ij}Z^0_{ab}+\frac2hM_{ij}Z^1_{ab}\right).
\]

`a^α·n=0`이므로 면내/두께 교차항은 없다. 각 길이 1인 `ζ` 구간에서
다음 `3×3` 행렬을 `(1,2,3)`와 `(3,4,5)`에 각각 더하면 `Z^0,Z^1`을 얻는다.

\[
Z^0_{\rm seg}=\frac1{30}\begin{bmatrix}4&2&-1\\2&16&2\\-1&2&4\end{bmatrix},\qquad
Z^1_{\rm seg}=\frac1{3}\begin{bmatrix}7&-8&1\\-8&16&-8\\1&-8&7\end{bmatrix}.
\]

따라서 `C=ρc(h/2)(M⊗Z0)`, `K=k[(h/2)(S⊗Z0)+(2/h)(M⊗Z1)]`이다.
이는 해당 node-major 순서와 일정 두께·물성·product metric에서만 적용한다.
요소/DOF source identity와 stable scatter를 보존하여 full `C,K,Q`를 조립한다.

## 5. 적분과 geometry 오차의 분리

Affine 평판(직사각형·평행사변형)의 일정 metric에서는 면내 `2×2` Gauss가 `M,S`에 정확하다.
두께 각 구간에서 `C`와 면내 전도는 4차 다항식이므로 **3점 Gauss**가 필요하다.
두께 전도는 2차이므로 2점으로 정확하지만 직접 검산은 통일된 3점을 사용한다.
이에 따라 affine의 직접 적분은 `2×2×(2 segments×3)=24`점으로 정확하다.
일반적인 평면 왜곡 Q4의 `S`는 역 Jacobian 때문에 다항식이 아니므로 이 exactness는 적용되지 않는다.

승인된 production 계산은 모든 요소의 면내 `3×3` Gauss와 위 exact `Z`의 분리 계산이다.
직접 계산하면 `3×3×2×3=54`점이며 같은 행렬이어야 한다. 3점 규칙은
`s=(-√(3/5),0,√(3/5))`, `w=(5/9,8/9,5/9)`이고 각 두께 구간에서 `dζ=ds/2`다.
표면의 warped/curved metric에는 유한 고정 차수의 exactness를 주장하지 않는다.
`3×3 → 5×5 → 7×7` 적분 수렴, mesh 수렴, thin-shell 모델 오차를 별도 확인한다.
3/5/7 차이는 진단용 evidence이며 3×3의 exactness 또는 전역 오차 bound가 아니다.
승인된 검증을 통과하지 못하면 조사·보고한다. 실패 후 지원 geometry 축소, 차수 변경,
adaptive quadrature 또는 tolerance 변경으로 통과시키지 않는다.

## 6. 절점 열입력

분포 edge 입력은 전부 연기되었다. 이전 초안의 edge 하중식·weight는 현재 kernel이나
검증 acceptance가 아니다. 입력 `Q`는 bottom/top의 signed inward nodal power만 포함한다.

Bottom/top 집중 열입력 `P [W]`는 지정 source node의 `a=1/5`에 그대로 더한다.
`h`, 면적 또는 적분 weight를 추가하지 않는다. 지정 온도와 같은 DOF의 입력도 full `Q`에
남겨 반력 계산에 반영한다. 중복 application의 가산/충돌 의미는 I/O 계약을 따른다.

## 7. Backward Euler, 지정 온도와 최초 증분

각 실제 증분 `δ_n>0`에서 backward Euler를 사용한다. `A_n=C/δ_n+K`이면

\[
C\frac{T^{n+1}-T^n}{\delta_n}+KT^{n+1}=Q+R^{n+1},\qquad
(A_n)_{ff}T_f^{n+1}=Q_f+\frac{C_{ff}}{\delta_n}T_f^n
+\frac{C_{fc}}{\delta_n}T_c^n-(A_n)_{fc}g.
\]

여기서 `f/c`는 free/constrained, `T_c^{n+1}=g`이다. 동등한 RHS 표현은
`Q_f−K_fc g+(C_ff/δ_n)T_f^n+(C_fc/δ_n)(T_c^n−g)`다.
**이전 constrained 온도 항은 첫 증분에서도 생략하지 않는다.** 이후 일정 `g`에서는
마지막 항이 사라진다. 잔차 tangent는 정확히 `A_n`; Newton 반복은 필요 없다.

`C>0`이므로 단열·무구속 모델도 transient `A_ff`는 양정치다. 정상 전도 `K`의 상수온도
null mode를 과도 해석의 singular error로 해석하지 않는다. 일반 partition 대수에서 모든
DOF가 constrained이면 `0×0` free solve를 건너뛰고 prescribed state와 full reaction을 계산한다.
이 경우는 generic partition fixture로만 검증한다. 이번 공개 입력은 bottom/top만 구속할 수
있으므로 내부 3개 station까지 모두 구속하는 열 CLI 모델은 허용되지 않는다.
고정 geometry·물성·DOF 집합과 동일한 `δ_n`이면 `A_ff` factorization을 재사용한다.
마지막 실제 `δ_n`이 바뀌면 다시 factorize한다. 근접한 계수라는 이유로 재사용하지 않는다.

**승인된 초기 상태:** 절점/nset 초기 온도를 5점에 복제한 전체 `T^0_raw`를 보존하며,
시간 적분의 `T^0=T^0_raw`로 사용한다. `t=0`에는 지정 온도로 덮어쓰지 않는다.
지정 온도는 첫 완료 증분 `t_1`부터 정확히 적용한다. 입력 초기값과 `g`가 다른 경우
그 차이는 첫 증분에 들어가므로 `C_fc(T_c^0−g)/δ_0`와 저장에너지가 모두 살아 있다.
이는 초기 불연속을 첫 유한 증분에서 처리하는 이산 정책으로, 실제 순간 impulse의
해상이나 `t=0+` 경계 만족을 주장하지 않는다. `t=0` 반력 열률은 시간 차분이 없어
미정의다. `t=0`에는 temperature만 저장하며 flux/reaction dataset을 만들지 않는다.
생략된 초기 온도는 0이다. Model-level BC를 먼저 적용하고 step-level BC가 같은 target을
override한다. 같은 stage 안에서 상충하는 duplicate 값은 오류다. 이 단계별 해석을
raw 초기 상태의 projection으로 바꾸지 않는다.

각 증분은 solve → full 온도 복구 → residual/flux/유한성·수지 확인 → candidate commit
순서로 처리한다. 실패 candidate가 이전 완료 state나 최종 HDF5를 대체하지 않는다.

### 7.1 재현 가능한 종료 시점과 증분 수

입력 binary64 `d=nominal Δt`, `τ=T_end`는 모두 유한한 양수다. `INC`는 양의 uint32,
기본값 100이며 허용 증분 수의 상한이다. `d≥τ`이면 한 증분, `δ_0=τ`, 저장 시각 `τ`다.
그 외에는 **반복 덧셈 없이** 정수 `j`의 곱으로 정상 시각을 만든다.

정수비 근방은 다음 연산 bound로 한 번만 판정한다. 먼저 비율이 `INC+1` 이내인지
확인하고 `τ`를 `[1/2,1)`에 놓는 같은 exact power-of-two로 `τ,d`를 재척도한다.
이 범위에서는 재척도된 `d`도 normal이며 아래 `p`의 overflow를 피할 수 있다.
아래 predicate는 재척도된 변수로 평가하며 결과와 시각·δ는 원래 단위를 유지한다.
`fl`은 binary64 nearest-even이고,
`rad(x)=max(x-nextafter(x,-∞), nextafter(x,+∞)-x)/2`는 해당 값의 보수적인
rounding-cell 반폭이다. 반폭이 subnormal로 표현되지 않으면 powers-of-two scale에서
계산한다. 여기서는 `q=fl(τ/d)`, `m=nearest_integer(q)`, `p=fl(m d)`를 쓴다.
`q`가 overflow하거나 `q>INC+1`이면 증분 수 초과다. 나머지 경우 `m≤2³²` 부근의
정수는 double에 정확히 표현된다. Tie의 정수 선택도 nearest-even으로 고정한다.

\[
B=d\,\operatorname{rad}(q)+\operatorname{rad}(p),\qquad
\widehat B=\operatorname{nextafter}(\operatorname{fl}(B),+\infty),\qquad
D=|\operatorname{fma}(-m,d,\tau)|,\qquad
\text{endpoint-close}\iff m\ge1\ \land\ D\le\widehat B.
\]

이는 임의 시간 epsilon이 아니다. 나눗셈 반올림의 dimensional uncertainty는
`d rad(q)`, 곱셈의 uncertainty는 `rad(p)`이므로 합이 두 연산에 대한 보수적 bound다.
특히 quotient가 정수 `m`으로 반올림되는 경우
`|τ-md|≤d rad(q)`이고, 정상 시각 곱 자체의 표시 불확실성은 `rad(p)`다.
Fused remainder로 곱의 선행 반올림을 제거하므로 `τ-fl(md)`의 이중 오차를 비교하지 않는다.
이 bound 안의 endpoint는 정상 마지막 증분과 같은 계수를 사용하는 명시적 정책이다.
실제 remainder가 이 범위 밖이면 작더라도 보존한다. 외부 입력 숫자의 임의 오차나 물리
시간 tolerance를 추가하지 않는다. 재척도된 `rad(q)`는 power-of-two이고 `d rad(q)`는
normal 범위의 exact binary scaling이다. 따라서 `fl(B)`는 덧셈 한 번의 반올림뿐이며
`nextafter` 한 번이 보수적 상계를 준다. FMA 한 번의 잔여 반올림까지 포함하면
accepted endpoint의 정확한 거리 bound는 `|τ-md|≤Bhat+rad(D)`다. 이렇게 **표현 가능한
threshold에 대한 `≤` 한 번**으로 판정하여 경계에서도 결과를 고정한다. FMA를 두 개의
일반 연산으로 분해하거나 compiler 재결합에 맡기지 않는다.

1. Endpoint-close이면 `N=m`. 아니면 `N=floor(q)+1=ceil(τ/d)`다. Quotient 반올림으로
   정수 경계를 넘을 수 있는 경우는 위 endpoint-close bound 안에 이미 포함되므로
   이 분기에서 floor는 안전하다. 먼저 nearest-integer branch를 판정하는 순서를 유지한다.
2. `N>INC`이면 solve 전에 실패한다. 정상 저장 시각은 `t_j=fl(j d)` (`j=1,…,N−1`),
   마지막은 **입력 그대로** `t_N=τ`다. 모든 시각이 유한하고 엄격히 증가하며
   `t_{N−1}<τ`인지 preflight에서 확인한다. `t_0=0`이다.
3. 정상 증분과 endpoint-close 마지막 증분의 연산 계수는 `δ=d`다. 실제로 짧아지는
   마지막 증분만 `δ=fl(τ-(N−1)d)`로 만든다. 이 remainder는
   `fma(-(N−1),d,τ)`처럼 곱-차를 한 번에 반올림한다. `δ>0`·유한성을 검사한다.
   저장된 `fl(jd)`끼리 다시 빼서 정상 `δ`를 만들지 않는다. 개념적인 격자 `j d`와
   저장 표현의 차이 및 endpoint-close 조정은 위 rounding bound로 제한된다.
4. 해당 증분의 `A`, RHS, full residual, reaction, `ΔE` 수지, BE 소산은 모두 **같은
   실제 `δ`**를 사용한다. Endpoint-close는 동일 계수를 재사용하며 실제 shortened
   마지막 `δ`가 `d`와 다를 때 refactorize한다.

예를 들어 `(τ,d)=(10,3)`은 `N=4`, `δ=(3,3,3,1)`, `(τ,d)=(2,3)`은 `N=1,δ=2`다.
`(0.9,0.3)` 같은 binary64 정수비 근방은 불필요한 미소 증분 없이 3회가 되며
세 번 모두 `δ=0.3`, 마지막 저장 시각은 `0.9`다. 표현오차만으로 재분해하지 않는다.
`nextafter`로 endpoint-close 경계를 양쪽에서 이동한 사례, 아주 작고 큰 시간 척도,
`INC=N`/`N−1`, uint32 overflow 방지를 독립 검증한다. 이 격자는 adaptive time stepping이 아니다.

## 8. 반력 열률, 에너지와 안정성

반력은 하중 제거 전 full residual에서 구한다.

\[
r^{n+1}=C(T^{n+1}-T^n)/\delta_n+KT^{n+1}-Q,\quad
R_c=r_c,\quad r_f\approx0,\qquad
U(T)=\mathbf1^TC(T-T_{\rm ref}\mathbf1).
\]

양의 `R_c [W]`는 지정 온도를 유지하려고 모델 **안으로** 공급하는 열률이다.
열원을 외부가 받는 부호는 반대다. 외부 reference의 부호는 별도 확인해야 한다.
`K1=0`과 대칭성에 의해 증분의 전역 수지는

\[
\Delta E=U^{n+1}-U^n=\delta_n\left(\mathbf1^TQ+\mathbf1_c^TR_c^{n+1}
+\mathbf1_f^Tr_f^{n+1}\right)
\]

이다. Physical balance defect `ΔE−δ_n(1ᵀQ+1_cᵀR_c)`와 `r_f`의 작음은 별도로
확인한다. 위 항등식에 `r_f`를 넣어 맞는다는 사실만으로 평형을 승인하지 않는다.
Full free residual과 energy는 내부 검증량이며 추가 공개 결과 quantity로 만들지 않는다.
정상상태에서는 총 입력과
반력 합이 0이다. 단열·무입력·무구속에서는 `U`가 보존되고 다음 이산 소산식이 성립한다.

\[
\tfrac12\big(\|T^{n+1}\|_C^2-\|T^n\|_C^2+\|T^{n+1}-T^n\|_C^2\big)
+\delta_n\|T^{n+1}\|_K^2=0.
\]

이 식의 이차량은 안정성 척도이며 실제 열에너지 `U`와 다르다.
Backward Euler는 이 의미에서 안정하지만 온도 positivity·maximum principle을 보장하지 않는다.
Consistent `C`, Q2 shape의 음의 구간과 일반 mesh의 `K`는 초기에 overshoot/undershoot를
만들 수 있다. 작은 `Δt`만으로 공간 이산의 단조성이 보장되지 않는다. Clamp, 임의 lumping,
인공 확산 또는 reference tolerance 완화로 이를 숨기지 않는다.

## 9. Heat flux 계산과 recovery identity

완료된 온도 `T_e^{n+1}`로 각 지정 평가점에서 다음을 직접 계산한다.

\[
q^{\rm global}=-kBT_e,\qquad
q^{\rm local}=\begin{bmatrix}e_1^T\\e_2^T\\n^T\end{bmatrix}q^{\rm global},\qquad
q_3=-\frac{2k}{h}\sum_{ia}N_iL_{a,\zeta}T_{ia}.
\]

Global 성분은 Cartesian `qx,qy,qz`, local은 `q1,q2,q3`이며 모두 열유속 `[W/m²]`다.
두께 적분 열류 `[W/m]`, 절점 입력 `[W]`, 반력 `[W]`와 구분한다.
예를 들어 `T=T_bottom+(ζ+1)(T_top−T_bottom)/2`면
`q3=−k(T_top−T_bottom)/h`; top이 뜨거우면 `q3<0`이다.
Bottom 외향 normal은 `−n`, top은 `+n`이므로 외향 성분은 각각 `−q3`, `+q3`다.

결정적인 필수 recovery 위치는 면내 4점
`(-g,-g),(g,-g),(g,g),(-g,g)`, `g=1/√3`와 두께 trace
`(-1,+),(-1/2,interior),(0,−),(0,+),(1/2,interior),(1,−)`의 곱인 **24점/요소**다.
두께 trace의 부호는 어느 구간 내부에서 접근하는지를 뜻한다. Integration 점과 recovery
점은 독립 계약이다. `ζ=0`의 온도는 같아도 `q3(0−),q3(0+)`는 서로 다를 수 있으며
이 둘을 평균하여 원자료를 잃지 않는다. 연속해의 재료 계면을 뜻하는 것도 아니다.

Row는 최소 `(instance, element label, in-plane point, thickness station, trace side, time)`와
component/좌표계로 구별한다. 공식 flux는 global Cartesian `q`이며 직교 frame
`R=[e1,e2,n]`을 보존하여 후처리에서 `q_local=Rᵀq`를 구한다.
요소 사이 edge의 flux도 요소별 한쪽 값을 유지한다. 서로 다른 `k`, normal, thickness에서
벡터를 무조건 같게 만들지 않는다. 접합부의 보존량은 약형으로 조립되는 열률이다.
`−kBT` raw flux의 경계 적분은 이산 반력과 일반적으로 정확히 같지 않으며 반력을 대신하지 않는다.
절점 외삽, nodal averaging, flux smoothing으로 원래 trace를 대체하지 않는다.

## 10. 검증 함의와 독립 검토 인계

다음은 실행 결과가 아닌 검증 설계다. 수치 기준은 독립 reviewer가 정하며
외부 비교는 [공통 TOLERANCE](../TOLERANCE.md)를 그대로 따른다.

| 검증 | 관찰할 내용 |
| --- | --- |
| Shape/행렬 | Nodal interpolation, partition of unity, `ΣB=0`, exact M/S/Z와 54점 직접 적분 일치, `C` SPD, `K` PSD/rank 19, symmetry, C·면내 K의 h scaling 및 두께 K의 1/h scaling |
| 평판 patch·회전 | Constant/면내 선형/두께 선형 온도에서 gradient·flux, 회전 시 T 불변과 q 회전, 법선 반전+station 재매핑 |
| 열입력 | Bottom/top 집중 입력 P 그대로, prescribed DOF와 겹친 P 보존과 반력 |
| 과도·초기 | 단열 균일온도 보존, 초기 불일치 첫 RHS, 공개 입력의 no-constrained, 에너지 수지·BE 소산; all-constrained는 generic partition fixture만 사용 |
| 면내 해석해 | 단열 직사각판 `T=T0+A cos(πx/L)exp[-α(π/L)²t]`, `α=k/(ρc)`; 면내 초기값은 5점 동일 |
| 두께 해석해 | 균일 초기온도 slab의 bottom 온도 step/top 단열에 대한 아래 series; 초기 jump와 첫 step 분리 |
| 재료·접속 | 두 물성/두께 strip의 직렬 전도는 아래 kernel/assembly oracle 범위로 구분; folded 등온/열률 수지와 junction 이상화의 별도 3D 한계 |
| 수렴·특이점 | 독립 semi-discrete eigenpair의 시간 1차 수렴, smooth affine mesh의 temperature L2 2차·gradient L2 1차, 집중 입력 근방 flux의 mesh 의존성과 곡률 오차 |
| Geometry·격자 | Whole-patch certificate, bow-tie/퇴화/uncertain metric, T-junction·비다양체·vertex-only fan·공유 edge 방향 오류, disconnected mesh, warped 3/5/7 적분 진단, 마지막 actual δ·endpoint-close·INC preflight |
| Reference | 기존 Step-1, `t=37`, NT11 bottom만 source identity로 비교; initial default·첫 증분·C 처리 차이 조사 |

단위 정사각형의 독립 exact surface oracle는 다음과 같다.

\[
M=\frac1{36}\begin{bmatrix}4&2&1&2\\2&4&2&1\\1&2&4&2\\2&1&2&4\end{bmatrix},\qquad
S=\frac1{6}\begin{bmatrix}4&-1&-2&-1\\-1&4&-1&-2\\-2&-1&4&-1\\-1&-2&-1&4\end{bmatrix}.
\]

`Z0·1=(1,4,2,4,1)/6`, `Z1·1=0`; 전체 thermal component마다 `K` null mode가
하나이며 모든 활성 DOF가 양의 체적 요소에 속하면 no-BC/disconnected 모델의
`C/δ+K`도 SPD다. 전체 K rank를 무조건 `N_dof−1`로 검사하지 않는다.
한쪽 trace fixture `(T1,…,T5)=(0,0,0,1,2)`는 `q3(0−)=0`, `q3(0+)=−4k/h`다.

최초 jump 독립 algebra fixture는 `C=[[2,1],[1,2]]`, `K=[[1,−1],[−1,1]]`,
`δ=1/2`, `Q=0`, raw `T0=(0,0)`, `g2=1`이다. 첫 두 증분의 free 온도는
`−1/5,1/25`, constrained reaction은 `24/5,36/25`이며 각각 `ΔE=δR_c`다.
이 기대값은 첫 `C_fc T_c_old/δ` 항을 검출한다. Generic all-constrained fixture와
함께 사용하되 공개 내부 station BC나 음수 clamp를 도입하지 않는다.

시간 수렴은 고정 mesh의 독립 generalized eigenpair `Kv=λCv`로
`T(t)=exp(−λt)v`와 BE의 `Π_j(1+λδ_j)⁻¹v`를 비교한다. 이 내부 초기 state를
공개 두께 비균일 초기 입력으로 확대하지 않는다. Temperature/gradient 공간 오차는
각각 `||T_h−T||_L2(V*)`, `||∇_*T_h−∇_*T||_L2(V*)`로 측정하고 시간 오차를 충분히
분리한다. 관찰 차수 `log(e_h/e_(h/2))/log(2)`의 수치 허용 구간은 reviewer가 고정한다.

Slab에서는 `x=z+h/2∈[0,h]`, bottom `Ts`, top 단열, 초기 `Ti`에 대해 `t>0`에서

\[
\theta=\frac{T-T_s}{T_i-T_s}
=\sum_{m=0}^{\infty}\frac4{(2m+1)\pi}
\sin\!\left(\frac{(2m+1)\pi x}{2h}\right)
\exp\!\left[-\alpha\left(\frac{(2m+1)\pi}{2h}\right)^2t\right],
\quad \alpha=\frac{k}{\rho c}.
\]

항별 미분으로 q를 구한다. 첫 생략 index가 `M`, `β=απ²t/(4h²)>0`이면 연속 항의
지수비는 최대 `r=exp[-8β(M+1)]<1`이다. 따라서 모든 `x`에서

\[
|\theta-\theta_M|\le\frac{4e^{-\beta(2M+1)^2}}{\pi(2M+1)(1-r)},\qquad
|q-q_M|\le\frac{2k|T_i-T_s|}{h}\frac{e^{-\beta(2M+1)^2}}{1-r}.
\]

`1-r`는 cancellation을 피하도록 `-expm1(-8β(M+1))`로 평가한다. 검증 시간·절단
`M`은 이 bound가 해당 oracle error budget보다 작도록 reviewer가 선택한다. 수치 적분한
slab energy는 temperature tail bound에 `ρc h·area·|Ti−Ts|`를 곱해 제한할 수 있다.

두 물성/두께 strip의 `R_th=Σ_e L_e/(k_e h_e w_e)` 정확해 검증은 양 끝 단면의 모든
station에 온도를 지정하는 kernel/assembly oracle이다. 내부 station 구속이 필요한 이
fixture는 공개 CLI acceptance 사례로 쓰지 않는다. 공개 입력으로 bottom/top만 구속한
strip은 두께 경계층을 포함한 독립 해 또는 동등한 조건의 oracle로 별도 검증해야 한다.

고정 5개 두께 DOF에서는 면내 refinement만으로 임의의 두께 경계층이 수렴하지 않는다.
Slab 해석해에 대한 두께 오차와 초기 `sqrt(αt)` 열침투 깊이를 확인해야 한다. 더 미세한
두께 모델은 독립 검증용 oracle이며 제품 입력의 온도점 수 확대를 뜻하지 않는다.
Reference CSV 한 시점은 다른 4개 온도점, flux, reaction 또는 전체 이력의 evidence가 아니다.
Consistent C와 본 적분·초기화 선택이 reference 정확도에 미치는 영향은 아직 측정되지 않았다.

Residual backward error는 `||r_f||/(||A_ff|| ||T_f||+||b_f||)`로, physical balance
defect는 `|ΔE−δ(1ᵀQ+1_cᵀR_c)|`를 저장·입력·반력 항 절댓값 합으로 정규화한다.
분모 0은 exact-zero case로 별도 처리한다. 행렬 오차도 해당 oracle norm과 연산 bound를
사용하며 단위 의존 `max(1,scale)`을 넣지 않는다. Geometry certificate·시간 predicate의
연산 bound와 convergence acceptance는 서로 다른 것이다.

Reviewer는 production kernel 착수 전에 roundoff/residual/energy 판정의 구체 bound,
mesh·시간열·series 절단·수렴차수 구간·warped quadrature 진단을 독립적으로 고정한다.
본 문서는 임의 tolerance 상수를 만들거나 공통 외부 5%를 내부 검증에 전용하지 않는다.
필수 NT11 비교는 `Step-1/t37/Part-1-1/source node/DOF11`의 729개 원본 row mapping과
[공통 TOLERANCE](../TOLERANCE.md)를 유지한다. 실패는 조사·상위 보고 대상이며
reference 보정·음수 clamp·방법 변경·tolerance 완화의 근거가 아니다.

물리·방법의 과거 미결정은 이번 승인 계약으로 대체되었다. 남은 독립 작업은 위 수치
acceptance 고정과 실제 구현 검증이다. Solver build/test, reference 수치 비교와 physics
검토는 본 정식화 작업에서 미실행이며, 스스로 수치 타당성이나 결과 correctness를 승인하지 않는다.