From e052b6e815698b00e032d2a2bde25d3145328bcb Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Tue, 29 Sep 2026 17:35:00 +0900 Subject: [PATCH] docs(thermal): specify approved formulation and robust geometry time grid --- .../formulation.md | 338 ++++++++++++++---- 1 file changed, 261 insertions(+), 77 deletions(-) diff --git a/docs/shell-transient-heat-transfer/formulation.md b/docs/shell-transient-heat-transfer/formulation.md index 2124e35..40a7f43 100644 --- a/docs/shell-transient-heat-transfer/formulation.md +++ b/docs/shell-transient-heat-transfer/formulation.md @@ -1,17 +1,19 @@ -# Shell 과도 열전달 정식화 제안 +# Shell 과도 열전달 정식화 -- 상태: **Proposed / 사용자 미승인 / 독립 numerical review 별도 보고서 참조**. 구현 실행 지시가 아니다. +- 상태: **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 내부 요소·적분·열용량·결과 생성과 동등하다는 주장이 아니다. +- 아래 수식은 승인된 FESA 모델이다. Abaqus 내부 알고리즘과 동등하다는 주장이 아니다. + 이전 Proposed 보고서의 검산은 제한된 과거 evidence이며 현재 solver 검증을 대신하지 않는다. ## 1. 적용 가정과 차원 고정된 3차원 midsurface의 4절점 Q4 요소에 두께 방향 온도점 5개를 둔다. -요소마다 두께 `h > 0`, 등방성 열전도율 `k > 0`, 밀도 `ρ > 0`, 비열 `c > 0`가 -상수이며 단면 전체는 같은 재료다. 서로 다른 요소에는 다른 값을 배정할 수 있다. -변형·대류·복사·내부 발열·적층·면적 분포 표면 열유속은 [PRD 제외 범위](PRD.md#제외-범위)를 따른다. +요소마다 유한한 두께 `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²]`이다. @@ -47,31 +49,110 @@ dV_* = J_s\frac h2\,d\xi\,d\eta\,d\zeta. 정확한 offset-volume의 접선 `aα+z n_,α`, 두께별 metric·면적 변화는 생략한다. 따라서 `h |κ1|`, `h |κ2|`가 작고 열장이 변하는 면내 길이에 비해 두께가 작은 영역을 -대상으로 한다. 적분 차수 증가로 이 모델링 오차가 없어지지는 않는다. 곡률의 정량 허용치, -왜곡 허용치와 validation margin은 아직 미정이다. 구조 shell의 두께 비율 기준을 전용하지 않는다. -곡률을 반영하는 별도 shell 열모델이 있다는 점은 [COMSOL의 곡면 열전달 예제](https://www.comsol.com/model/heat-transfer-in-curved-layered-shells-104241)로 확인할 수 있다. +대상으로 하는 물리 근사다. 이는 승인되지 않은 curvature/thickness cutoff를 입력에 추가하는 +조건이 아니다. 평판·warped patch·곡면 mesh·fold를 지원하며 적분 차수를 높여도 생략한 +offset-volume 모델링 오차는 없어지지 않는다. Geometry의 필수 조건은 유한 좌표, 서로 다른 4절점, `J_s>0`, 비특이 metric, 자기 교차 없는 유효 patch와 일관된 orientable connectivity이다. 적분점 검사만으로 전체 patch 유효성을 -증명하지 않는다. 검사 위치·scale-aware degeneracy 기준·비다양체 허용 범위는 I/O와 -numerical review가 확정해야 하며 작은 양수 Jacobian을 clamp하지 않는다. +증명하지 않는다. 아래 whole-patch certificate를 사용하며 작은 양수 Jacobian을 clamp하지 않는다. -**접속 제안:** 인접 요소는 같은 source node의 같은 `ζ` 온도 DOF를 공유한다. +**승인된 접속:** 인접 요소는 같은 source node의 같은 `ζ` 온도 DOF를 공유한다. 매끈한 면에서 이는 온도 연속성과 약형 열률 평형을 준다. 접힌 판에서는 요소별 normal을 유지하고 공유 edge의 `ζ`별 온도만 같게 묶는 이상적 무저항·무용량 접합이다. 서로 다른 두께에서는 같은 normalized `ζ`를 연결하므로 실제 `z` 좌표가 같다는 뜻이 아니다. 실제 step 단차의 노출면, 유한 접합부 부피, 국소 spreading resistance는 해상하지 않는다. -이 모델은 사용자 결정이 필요한 물리 이상화다. 동의 없이 실제 3D 접속을 지원한다고 표시하지 않는다. +이는 승인된 이상화이며 상세 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의 물리 위치를 -혼동하지 않는다. 이 외부 위치 의미는 [Abaqus Heat Transfer Shell Elements](https://docs.software.vt.edu/abaqusv2025/English/SIMACAEELMRefMap/simaelm-c-shellelem.htm)의 공개 규약을 따른다. +혼동하지 않는다. 외부 위치 의미는 승인된 DS4 입력 계약을 따른다. 면내 shape와 미분은 각 절점의 부호 `(ξ_i,η_i)`를 사용한다. @@ -81,7 +162,7 @@ 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 구간**으로 제안한다. +두께는 `I−=[-1,0]`, `I+=[0,1]`의 **두 C0 quadratic 구간**을 사용한다. 단일 quartic보다 각 구간의 저차 보간·국소 미분을 명료하게 유지하며, 중간 온도 1개를 공유한다. 구간 좌표는 `s=2ζ+1` 또는 `s=2ζ−1`; 두 경우 `s∈[-1,1]`, `ds/dζ=2`이다. @@ -113,7 +194,7 @@ Fourier 법칙 `q=−k∇_*T`, 저장식 `ρc Tdot + div_* q=0`를 사용한다. \[ \int_{V_*}\rho c\,v\dot T\,dV_*+ \int_{V_*}k\nabla_*v\cdot\nabla_*T\,dV_* -=\int_{\Gamma_{\rm edge,*}}v q_{\rm in}\,dA_*+\sum_j v_jP_j. +=\sum_j v_jP_j. \] 지정 온도 외의 하중 없는 경계는 자연 단열이다. 집중 입력 `P_j`는 승인된 bottom/top @@ -127,7 +208,7 @@ Fourier 법칙 `q=−k∇_*T`, 저장식 `ρc Tdot + div_* q=0`를 사용한다. `C_e [J/K]`는 consistent 열용량, `K_e [W/K]`는 열전도 행렬이다. `C_e` 대칭 양정치, `K_e` 대칭 양반정치와 유일한 요소 상수온도 null mode를 기대한다. -Reduced integration, hourglass stiffness, capacity lumping은 이 제안에 넣지 않는다. +Reduced integration, hourglass stiffness, capacity lumping, 인공 확산은 사용하지 않는다. 명시적인 component식과 분리 계산에 필요한 작은 행렬은 다음과 같다. @@ -151,93 +232,136 @@ 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,F`를 조립한다. +요소/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점을 권장한다. +두께 전도는 2차이므로 2점으로 정확하지만 직접 검산은 통일된 3점을 사용한다. 이에 따라 affine의 직접 적분은 `2×2×(2 segments×3)=24`점으로 정확하다. 일반적인 평면 왜곡 Q4의 `S`는 역 Jacobian 때문에 다항식이 아니므로 이 exactness는 적용되지 않는다. -제안 기본은 모든 요소의 면내 `3×3` Gauss와 위 exact `Z`의 분리 계산이다. +승인된 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 모델 오차를 별도 확인한다. -지원 geometry 범위에서 기본 차수가 불충분하면 검토를 통해 차수/범위를 고정하고 구현한다. -실행마다 허용오차에 맞춰 임의 adaptive quadrature를 선택하는 계약은 아니다. +3/5/7 차이는 진단용 evidence이며 3×3의 exactness 또는 전역 오차 bound가 아니다. +승인된 검증을 통과하지 못하면 조사·보고한다. 실패 후 지원 geometry 축소, 차수 변경, +adaptive quadrature 또는 tolerance 변경으로 통과시키지 않는다. -## 6. 가장자리·절점 열입력 +## 6. 절점 열입력 -Edge 순서는 `1:(1,2), 2:(2,3), 3:(3,4), 4:(4,1)`로 제안한다. -면내 자연좌표 edge curve를 `r_edge(u)`, `u∈[-1,1]`, `J_l=|dr_edge/du|`로 둔다. -두께 전체 균일 입력 `q_in [W/m²]`의 일관 하중은 다음과 같다. - -\[ -(F_e^{\rm edge})_{ia}=\int_{-1}^1\int_{-1}^1 -N_i(\xi(u),\eta(u))L_a(\zeta)q_{\rm in}J_l\frac h2\,d\zeta du. -\] - -Q4의 직선 edge, 일정 입력·두께에는 총 입력 `q_in h L_edge`가 정확히 보존된다. -각 endpoint의 5성분은 `q_in h L_edge (1,4,2,4,1)/24`이며 나머지 두 절점은 0이다. -이는 `∫L dζ=(1,4,2,4,1)/6`와 edge당 `∫N_i dl=L_edge/2`에서 나온다. -Edge 2점, 각 두께 구간 2점으로 이 하중은 정확하다. 내부 온도점에도 일관 하중이 생기는 -것은 허용된 edge 입력의 결과이며, 내부 온도점 직접 집중 입력을 지원한다는 뜻은 아니다. -외부 syntax가 선열률 `λ [W/m]`를 쓰기로 결정되면 `q_in=λ/h`로 해석해야 한다. -두 단위를 혼용하지 않으며 외부 입력 단위·keyword는 I/O owner가 확정한다. +분포 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 `F`에 +`h`, 면적 또는 적분 weight를 추가하지 않는다. 지정 온도와 같은 DOF의 입력도 full `Q`에 남겨 반력 계산에 반영한다. 중복 application의 가산/충돌 의미는 I/O 계약을 따른다. ## 7. Backward Euler, 지정 온도와 최초 증분 -고정 `Δt>0`에서 backward Euler를 제안한다. `t_n=nΔt`, `A=C/Δt+K`이면 +각 실제 증분 `δ_n>0`에서 backward Euler를 사용한다. `A_n=C/δ_n+K`이면 \[ -C\frac{T^{n+1}-T^n}{\Delta t}+KT^{n+1}=F^{n+1}+R^{n+1},\qquad -A_{ff}T_f^{n+1}=F_f^{n+1}+\frac{C_{ff}}{\Delta t}T_f^n -+\frac{C_{fc}}{\Delta t}T_c^n-A_{fc}g^{n+1}. +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^{n+1}`이다. 동등한 RHS 표현은 -`F_f−K_fc g^{n+1}+(C_ff/Δt)T_f^n+(C_fc/Δt)(T_c^n−g^{n+1})`다. +여기서 `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`; Newton 반복은 필요 없다. -Backward Euler의 mass/diffusion 행렬 구성은 [MFEM Example 16 원문](https://raw.githubusercontent.com/mfem/mfem/master/examples/ex16.cpp)에서도 확인할 수 있으며, 위 제약·초기화 정책은 별도의 FESA 제안이다. +마지막 항이 사라진다. 잔차 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 집합·Δt에서는 `A_ff`를 최초 한 번 factorize할 수 있다. -전체 시간이 `Δt`의 정수 배가 아니면 입력을 거부하는 방안을 권장한다. 짧은 마지막 증분은 -고정 `Δt`와 factorization 재사용 의미를 바꾸므로 자동 도입하지 않는다. 시간 비교의 수치 -판정과 종료 규칙은 사용자 검토와 I/O 계약에서 확정한다. +고정 geometry·물성·DOF 집합과 동일한 `δ_n`이면 `A_ff` factorization을 재사용한다. +마지막 실제 `δ_n`이 바뀌면 다시 factorize한다. 근접한 계수라는 이유로 재사용하지 않는다. -**초기 불일치 제안:** 입력 초기 온도를 5점에 복제한 전체 `T^0_raw`를 보존하며, +**승인된 초기 상태:** 절점/nset 초기 온도를 5점에 복제한 전체 `T^0_raw`를 보존하며, 시간 적분의 `T^0=T^0_raw`로 사용한다. `t=0`에는 지정 온도로 덮어쓰지 않는다. 지정 온도는 첫 완료 증분 `t_1`부터 정확히 적용한다. 입력 초기값과 `g`가 다른 경우 -그 차이는 첫 증분에 들어가므로 `C_fc(T_c^0−g)/Δt`와 저장에너지가 모두 살아 있다. +그 차이는 첫 증분에 들어가므로 `C_fc(T_c^0−g)/δ_0`와 저장에너지가 모두 살아 있다. 이는 초기 불연속을 첫 유한 증분에서 처리하는 이산 정책으로, 실제 순간 impulse의 해상이나 `t=0+` 경계 만족을 주장하지 않는다. `t=0` 반력 열률은 시간 차분이 없어 -미정의다. 0으로 만들어 ‘계산된 반력’으로 저장하지 않는다. -대안인 `t=0` 즉시 projection은 별도 에너지 jump `1ᵀC(T^0_projected−T^0_raw)`와 -그 보존/보고 정책이 필요하다. 이 선택은 결과를 바꾸므로 사용자 결정 후 하나만 구현한다. -생략된 초기 온도의 값은 이 문서에서 결정하지 않는다. Reference는 이 default에 민감하다. +미정의다. `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 t+KT^{n+1}-F^{n+1},\quad +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). \] @@ -247,16 +371,19 @@ U(T)=\mathbf1^TC(T-T_{\rm ref}\mathbf1). `K1=0`과 대칭성에 의해 증분의 전역 수지는 \[ -U^{n+1}-U^n=\Delta t\left(\mathbf1^TF^{n+1}+\mathbf1_c^TR_c^{n+1} +\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) \] -이다. Free residual을 버리지 않고 수지 defect와 함께 확인한다. 정상상태에서는 총 입력과 +이다. 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 t\|T^{n+1}\|_K^2=0. ++\delta_n\|T^{n+1}\|_K^2=0. \] 이 식의 이차량은 안정성 척도이며 실제 열에너지 `U`와 다르다. @@ -281,7 +408,7 @@ Global 성분은 Cartesian `qx,qy,qz`, local은 `q1,q2,q3`이며 모두 열유 `q3=−k(T_top−T_bottom)/h`; top이 뜨거우면 `q3<0`이다. Bottom 외향 normal은 `−n`, top은 `+n`이므로 외향 성분은 각각 `−q3`, `+q3`다. -결정적인 필수 recovery 위치 제안은 면내 4점 +결정적인 필수 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 @@ -289,29 +416,77 @@ Bottom 외향 normal은 `−n`, top은 `+n`이므로 외향 성분은 각각 ` 이 둘을 평균하여 원자료를 잃지 않는다. 연속해의 재료 계면을 뜻하는 것도 아니다. Row는 최소 `(instance, element label, in-plane point, thickness station, trace side, time)`와 -component/좌표계로 구별해야 한다. Global/local 중 영속 저장 형식은 I/O owner가 결정한다. +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은 이 제안의 필수 recovery가 아니다. +절점 외삽, nodal averaging, flux smoothing으로 원래 trace를 대체하지 않는다. -## 10. 검증 함의와 미결정 사항 +## 10. 검증 함의와 독립 검토 인계 다음은 실행 결과가 아닌 검증 설계다. 수치 기준은 독립 reviewer가 정하며 외부 비교는 [공통 TOLERANCE](../TOLERANCE.md)를 그대로 따른다. | 검증 | 관찰할 내용 | | --- | --- | -| Shape/행렬 | Nodal interpolation, partition of unity, `ΣB=0`, explicit Z와 직접 적분 일치, `C` SPD, `K` PSD/rank 19, symmetry | +| 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 재매핑 | -| 열입력 | Edge 총률 `q_in hL`·5점 weight, bottom/top 집중 입력, prescribed DOF와 겹친 입력의 반력 | +| 열입력 | 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의 두 표면 온도 step에 대한 1D Fourier-series 해와 비교; 초기 jump와 첫 step 분리 | +| 두께 해석해 | 균일 초기온도 slab의 bottom 온도 step/top 단열에 대한 아래 series; 초기 jump와 첫 step 분리 | | 재료·접속 | 두 물성/두께 strip의 직렬 전도는 아래 kernel/assembly oracle 범위로 구분; folded 등온/열률 수지와 junction 이상화의 별도 3D 한계 | -| 수렴·특이점 | `Δt` 1차 수렴, 면내 mesh/적분 차수 수렴, 집중 입력 근방 flux의 mesh 의존성과 곡률 오차 | +| 수렴·특이점 | 독립 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만 구속한 @@ -323,10 +498,19 @@ Slab 해석해에 대한 두께 오차와 초기 `sqrt(αt)` 열침투 깊이를 Reference CSV 한 시점은 다른 4개 온도점, flux, reaction 또는 전체 이력의 evidence가 아니다. Consistent C와 본 적분·초기화 선택이 reference 정확도에 미치는 영향은 아직 측정되지 않았다. -미결정 owner는 메인/사용자(물리·초기 정책), I/O owner(공개 의미), numerical reviewer(수치 허용 범위)다. -사용자 검토에는 thin-shell product geometry와 접속 이상화, C/BE, 최초 불일치·종료 시점, -24개 flux trace 제안을 포함한다. Geometry 판정 기준과 기본 적분의 지원 범위, 생략 초기값, -edge 외부 단위, t=0 결과의 validity 표현은 구현 전에 확정해야 한다. -현재 정식화 담당자의 확인은 문헌·요구사항 대조와 문서 정적 검토뿐이다. Solver build/test, -reference 수치 비교와 physics 검토는 미실행이다. 독립 numerical review의 결과는 별도 -보고서를 참조하며 본 문서는 스스로 수치 타당성이나 결과 correctness를 승인하지 않는다. +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를 승인하지 않는다.