FESADev/docs/formulations/3d-isoparametric-euler-beam-formulation.md

3D Isoparametric Euler–Bernoulli Beam Formulation

Metadata

• feature_id: 3d-isoparametric-euler-beam
• source_requirement: docs/requirements/linear-static-3d-euler-beam.md
• source_research: docs/research/linear-static-3d-euler-beam-research.md
• approved_design: docs/superpowers/specs/2026-08-08-linear-static-3d-euler-beam-design.md
• status: ready-for-numerical-review
• owner_agent: formulation-agent
• date: 2026-08-09

1. 목적과 정식화 요약

이 문서는 선형 정적 구조해석 파이프라인을 검증하기 위한 2절점 3차원 Euler–Bernoulli beam 요소의 수학 계약을 정의한다. 요소는 축 변형, Saint-Venant 비틀림, 두 주축 방향 굽힘을 포함하며 절점마다 3개의 병진 자유도와 3개의 회전 자유도를 가진다.

요소 계산의 핵심 흐름은 다음과 같다.

\boxed{
\mathbf X_1,\mathbf X_2,\mathbf a
\rightarrow L,\mathbf R,\mathbf T
\rightarrow \mathbf B(\xi),\mathbf D
\rightarrow \mathbf K_e^l
\rightarrow \mathbf K_e^g
\rightarrow \mathbf K\mathbf d=\mathbf F
\rightarrow \boldsymbol\varepsilon,\mathbf s,\mathbf f_e
}

행렬식으로 쓰면

\begin{aligned}
\boldsymbol\varepsilon(\xi)&=\mathbf B(\xi)\mathbf d_e^l,\\
\mathbf s(\xi)&=\mathbf D\boldsymbol\varepsilon(\xi),\\
\mathbf K_e^l&=\int_0^L\mathbf B^T\mathbf D\mathbf B\,dx,\\
\mathbf K_e^g&=\mathbf T^T\mathbf K_e^l\mathbf T,\\
\mathbf r_e^g&=\mathbf K_e^g\mathbf d_e^g-\mathbf f_e^g.
\end{aligned}

성분형은

\varepsilon_i=\sum_{a=1}^{12}B_{ia}d_a^l,\qquad
s_i=\sum_{j=1}^{4}D_{ij}\varepsilon_j,

K_{ab}^{l}=\int_0^L\sum_{i=1}^{4}\sum_{j=1}^{4}
B_{ia}D_{ij}B_{jb}\,dx,\qquad
r_a^g=\sum_{b=1}^{12}K_{ab}^{g}d_b^g-f_a^g.

Beam and Frame Finite Elements는 Euler–Bernoulli beam의 횡변위와 회전 자유도에 cubic Hermite 보간이 필요하고, 공간 frame 요소가 축·굽힘·비틀림과 좌표 변환을 결합한다는 이론적 근거를 제공한다. Isoparametric Finite Elements는 자연좌표, Jacobian, 미분 변환 및 수치적분의 공통 틀을 제공한다.

2. Scope and Assumptions

2.1 포함 범위

• 2절점 직선 prismatic 공간 beam 요소
• Abaqus TYPE=B33에만 연결되는 shear-free cubic Euler–Bernoulli 의미
• small displacement 및 small rotation
• 선형 정적 해석
• 균질 선형 탄성 재료
• 축 변형, 두 방향 Euler–Bernoulli 굽힘, Saint-Venant 비틀림
• 절점당 자유도 [u,v,w,\theta_x,\theta_y,\theta_z]
• 단면 도심을 지나는 국부 $x$축
• 단면 주축과 국부 $y,z$축의 일치: I_{yz}=0
• 보존력과 일관 등가절점하중
• formulation 검증용 상수 국부 line-load 등가절점벡터; Abaqus *DLOAD 입력은 제외
• 자연좌표 $\xi\in[-1,1]$을 사용하는 1차원 등매개 기하 mapping

2.2 제외 범위

• Timoshenko 전단 변형과 shear correction factor
• Abaqus TYPE=B31의 Euler–Bernoulli 요소 재해석
• 기하비선형, 재료비선형, P–\Delta, geometric stiffness
• 단면 warping, bimoment 및 warping 자유도
• shear-center/centroid offset, rigid offset 및 end release
• 곡선 beam, taper 또는 요소 내부에서 변하는 E,G,A,I_y,I_z,J
• 비주축 단면의 I_{yz}\ne0 굽힘 결합
• 질량, 감쇠, 동적 해석
• 열변형, 초기변형률, 초기응력
• follower load

2.3 “Isoparametric”의 정확한 의미

이 요소에서 두 절점의 위치와 축변위 및 비틀림 회전은 같은 2절점 Lagrange 함수로 보간한다. Euler–Bernoulli 굽힘은 횡변위와 그 기울기를 동시에 절점에서 보간해야 하므로 cubic Hermite 함수를 사용한다. 따라서 이 문서에서 “isoparametric”은 자연좌표 기반의 기하 mapping과 적분 틀을 뜻하며, 여섯 운동학 변수 모두에 동일한 Lagrange 함수를 사용한다는 뜻은 아니다.

2.4 단위

일반 입력에는 특정 단위계를 고정하지 않고 일관 단위계를 사용한다. 승인된 reference/cantilever beam/ 비교만 외부 계약에 따라 SI를 사용하며, 이 사실은 일반 .inp 입력의 단위계를 SI로 추론한다는 뜻이 아니다.

양	차원
x,L,u,v,w	length
\theta_x,\theta_y,\theta_z	dimensionless (radian)
E,G	force / length$^2$
A	length$^2$
I_y,I_z,J	length$^4$
N	force
T,M_y,M_z	force \times length

3. 좌표계와 부호 규약

3.1 국부 직교좌표계

전역 절점 위치를 $\mathbf X_1,\mathbf X_2\in\mathbb R^3$라 하고 Abaqus first section axis n1으로 제공되는 guide vector를 $\mathbf a=\mathbf n_1$이라 한다. 요소 tangent $\mathbf t$를 국부 $x$축으로, 투영·정규화한 n1을 국부 $y$축으로, \mathbf t\times\mathbf n_1 방향을 국부 $z$축으로 둔다.

\begin{aligned}
\Delta\mathbf X&=\mathbf X_2-\mathbf X_1,\\
L&=\|\Delta\mathbf X\|,\\
\mathbf e_x&=\frac{\Delta\mathbf X}{L},\\
\widetilde{\mathbf e}_y&=\mathbf a-(\mathbf a\cdot\mathbf e_x)\mathbf e_x,\\
\mathbf e_y&=\frac{\widetilde{\mathbf e}_y}
{\|\widetilde{\mathbf e}_y\|},\\
\mathbf e_z&=\mathbf e_x\times\mathbf e_y.
\end{aligned}

성분형은

\Delta X_I=X_{2I}-X_{1I},\qquad
L=\sqrt{\sum_{I=1}^{3}(\Delta X_I)^2},\qquad
e_{xI}=\frac{\Delta X_I}{L},

\widetilde e_{yI}=a_I-\left(\sum_{J=1}^{3}a_Je_{xJ}\right)e_{xI},
\qquad
e_{yI}=\frac{\widetilde e_{yI}}
{\sqrt{\sum_{J=1}^{3}\widetilde e_{yJ}^2}},

e_{zI}=\sum_{J=1}^{3}\sum_{K=1}^{3}\epsilon_{IJK}e_{xJ}e_{yK}.

좌표 크기와 guide vector 크기를 반영한 승인된 validation은 다음과 같다.

L>10^{-12}\max(1,\|\mathbf X_1\|,\|\mathbf X_2\|),\qquad
\|\widetilde{\mathbf e}_y\|>
10^{-12}\max(1,\|\mathbf a\|),\qquad
\mathbf e_i\cdot\mathbf e_j=\delta_{ij},\qquad
\mathbf e_x\times\mathbf e_y=\mathbf e_z.

zero guide vector, guide vector가 요소 축과 평행한 경우, 또는 zero/near-zero length는 국부 $y,z$축과 $I_y,I_z$의 의미를 정할 수 없으므로 model validation 오류로 처리한다. 임의의 fallback 축을 조용히 선택하지 않는다.

3.2 변위와 회전 부호

• u,v,w: 각각 국부 +x,+y,+z 방향 병진변위
• \theta_x,\theta_y,\theta_z: 각각 국부축에 대한 오른손 법칙의 양의 회전
• Euler–Bernoulli 무전단 조건:

\boxed{\theta_z=\frac{dv}{dx},\qquad
\theta_y=-\frac{dw}{dx}}

두 식의 부호가 다르다는 점이 $w$–\theta_y 굽힘 강성행렬의 6EI/L^2 항 부호를 결정한다.

3.3 단면력 부호

단면 일반화 변형률과 단면력을 다음 순서로 둔다.

\boldsymbol\varepsilon=
\begin{bmatrix}\varepsilon_0&\kappa_x&\kappa_y&\kappa_z\end{bmatrix}^T,
\qquad
\mathbf s=
\begin{bmatrix}N&T&M_y&M_z\end{bmatrix}^T.

$N$은 인장일 때 양수이고 $T,M_y,M_z$는 해당 국부축 오른손 법칙을 따른다.

4. Primary Variables and DOFs

4.1 연속장

약형의 독립장은 중심선 축변위 u(x), 비틀림 회전 \theta_x(x), 두 횡변위 $v(x),w(x)$이다. 단면 굽힘 회전은 횡변위 기울기에서 유도한다.

u,\theta_x\in H^1(0,L),\qquad v,w\in H^2(0,L).

4.2 절점 자유도 벡터

국부 요소 자유도 순서는 다음으로 고정한다.

\boxed{
\mathbf d_e^l=
\begin{bmatrix}
u_1&v_1&w_1&\theta_{x1}&\theta_{y1}&\theta_{z1}&
u_2&v_2&w_2&\theta_{x2}&\theta_{y2}&\theta_{z2}
\end{bmatrix}^T
}

성분 번호는 다음과 같다.

a	1	2	3	4	5	6	7	8	9	10	11	12
d_a^l	u_1	v_1	w_1	\theta_{x1}	\theta_{y1}	\theta_{z1}	u_2	v_2	w_2	\theta_{x2}	\theta_{y2}	\theta_{z2}

전역 nodal component 순서는 각 node에서 [UX, UY, UZ, URX, URY, URZ]이고, $\mathbf d_e^l=\mathbf T\mathbf d_e^g$는 이 node/component block 순서를 유지한 채 국부 $[u,v,w,\theta_x,\theta_y,\theta_z]$로 회전한다.

요소는 equation id를 소유하지 않는다. 전역 equation numbering과 constrained/free mapping은 프로젝트 아키텍처에 따라 DofManager가 소유해야 한다.

5. Strong Form and Boundary Conditions

5.1 단면력–변위 관계

상수 단면과 선형 탄성에서

\begin{aligned}
N&=EA\,u',\\
T&=GJ\,\theta_x',\\
M_y&=EI_y\,\theta_y'=-EI_y\,w'',\\
M_z&=EI_z\,\theta_z'=EI_z\,v''.
\end{aligned}

여기서 $(\cdot)'=d(\cdot)/dx$이다.

5.2 미분방정식

국부 분포하중을 p_x,p_y,p_z 및 단위길이당 분포 비틀림 모멘트를 $m_x$라 하면 요소 내부 평형식은

\boxed{
\begin{aligned}
-N'&=p_x,\\
-T'&=m_x,\\
M_z''&=p_y,\\
-M_y''&=p_z.
\end{aligned}}

변위로 직접 쓰면

\boxed{
\begin{aligned}
-(EAu')'&=p_x,\\
-(GJ\theta_x')'&=m_x,\\
(EI_zv'')''&=p_y,\\
(EI_yw'')''&=p_z.
\end{aligned}}

상수 물성이고 요소 내부 분포하중이 없으면

u''=0,\qquad \theta_x''=0,\qquad v''''=0,\qquad w''''=0.

따라서 선형 축/비틀림 보간과 cubic Hermite 굽힘 보간이 각각 homogeneous strong form을 정확히 만족한다.

5.3 경계조건과 자연경계 부호

경계 외향 법선의 1차원 부호를 n=-1 at x=0, n=+1 at $x=L$로 둔다. 각 끝에서는 다음 essential/natural 쌍 중 하나를 지정한다.

장	essential	natural
축	u=\bar u	F_x=nN
비틀림	\theta_x=\bar\theta_x	M_x=nT
$x$–y 평면 굽힘	v=\bar v	F_y=-nM_z'
$x$–y 평면 회전	\theta_z=\bar\theta_z	M_z^{\mathrm{end}}=nM_z
$x$–z 평면 굽힘	w=\bar w	F_z=nM_y'
$x$–z 평면 회전	\theta_y=\bar\theta_y	M_y^{\mathrm{end}}=nM_y

굽힘 방정식은 4차이므로 각 굽힘 평면에서 변위/회전 또는 전단력/모멘트의 두 경계조건이 필요하다.

6. Weak or Variational Form

허용 시험함수 $\delta u,\delta\theta_x,\delta v,\delta w$는 essential boundary에서 0이다. 두 번의 부분적분을 적용하면 다음 약형을 얻는다.

\boxed{
\begin{aligned}
&\int_0^L
\left[
EA\,\delta u' u'
+GJ\,\delta\theta_x'\theta_x'
+EI_y\,\delta w''w''
+EI_z\,\delta v''v''
\right]dx\\
&=\int_0^L
\left[
\delta u\,p_x
+\delta\theta_x\,m_x
+\delta v\,p_y
+\delta w\,p_z
\right]dx
+\delta\mathbf d_{\Gamma}^T\bar{\mathbf f}_{\Gamma}.
\end{aligned}}

경계 가상일의 성분형은

\delta\mathbf d_{\Gamma}^T\bar{\mathbf f}_{\Gamma}
=\sum_{x\in\{0,L\}}
\left(
\delta u F_x+\delta v F_y+\delta w F_z
+\delta\theta_xM_x+\delta\theta_yM_y^{\mathrm{end}}
+\delta\theta_zM_z^{\mathrm{end}}
\right).

일반화 변형률과 단면력을 사용하면 내부 가상일은 더 간단히

\delta W_{\mathrm{int}}^e
=\int_0^L\delta\boldsymbol\varepsilon^T\mathbf s\,dx
=\int_0^L\delta\boldsymbol\varepsilon^T
\mathbf D\boldsymbol\varepsilon\,dx

이고, 성분형은

\delta W_{\mathrm{int}}^e
=\int_0^L
\left(
N\,\delta\varepsilon_0
+T\,\delta\kappa_x
+M_y\,\delta\kappa_y
+M_z\,\delta\kappa_z
\right)dx.

선형 탄성에서는 동일한 식을 전체 잠재에너지의 정류조건으로 쓸 수 있다.

\Pi_e(\mathbf d_e^l)
=\frac12(\mathbf d_e^l)^T\mathbf K_e^l\mathbf d_e^l
-(\mathbf d_e^l)^T\mathbf f_e^l,
\qquad
\frac{\partial\Pi_e}{\partial\mathbf d_e^l}=\mathbf0.

7. Mapping and Natural Coordinates

7.1 2절점 등매개 기하 mapping

자연좌표 $\xi\in[-1,1]$에서

N_1^L(\xi)=\frac{1-\xi}{2},\qquad
N_2^L(\xi)=\frac{1+\xi}{2}.

기하 위치는

\mathbf X(\xi)=N_1^L(\xi)\mathbf X_1+N_2^L(\xi)\mathbf X_2.

국부 축좌표 $x\in[0,L]$는

x(\xi)=\frac{L}{2}(1+\xi),\qquad
r=\frac{x}{L}=\frac{1+\xi}{2}.

Jacobian과 미분 변환은

J_x=\frac{dx}{d\xi}=\frac{L}{2}>0,
\qquad
dx=J_xd\xi,

\frac{d}{dx}=\frac{1}{J_x}\frac{d}{d\xi}
=\frac{2}{L}\frac{d}{d\xi},\qquad
\frac{d^2}{dx^2}=\frac{4}{L^2}\frac{d^2}{d\xi^2}.

직선 2절점 요소에서는 $J_x$가 상수이다. $J_x\le J_{\mathrm{tol}}$이면 zero 또는 near-zero length 요소이므로 계산을 중단한다.

8. Discretization

8.1 축변위와 비틀림 보간

u(\xi)=N_1^Lu_1+N_2^Lu_2,
\qquad
\theta_x(\xi)=N_1^L\theta_{x1}+N_2^L\theta_{x2}.

성분형 보간 조건은

N_i^L(\xi_j)=\delta_{ij},\qquad
N_1^L+N_2^L=1.

8.2 Cubic Hermite 함수

$r=x/L$을 사용하면

\begin{aligned}
H_1(r)&=1-3r^2+2r^3,\\
H_2(r)&=L(r-2r^2+r^3),\\
H_3(r)&=3r^2-2r^3,\\
H_4(r)&=L(-r^2+r^3).
\end{aligned}

자연좌표 성분형은

\begin{aligned}
H_1(\xi)&=\frac14(2-3\xi+\xi^3),\\
H_2(\xi)&=\frac{L}{8}(1-\xi-\xi^2+\xi^3),\\
H_3(\xi)&=\frac14(2+3\xi-\xi^3),\\
H_4(\xi)&=\frac{L}{8}(-1-\xi+\xi^2+\xi^3).
\end{aligned}

횡변위 보간은 회전 부호를 반영하여

\boxed{
\begin{aligned}
v(x)&=H_1v_1+H_2\theta_{z1}+H_3v_2+H_4\theta_{z2},\\
w(x)&=H_1w_1-H_2\theta_{y1}+H_3w_2-H_4\theta_{y2}.
\end{aligned}}

1차 미분은

\begin{aligned}
H_1'&=\frac{6(r^2-r)}{L},&
H_2'&=1-4r+3r^2,\\
H_3'&=\frac{6(r-r^2)}{L},&
H_4'&=-2r+3r^2.
\end{aligned}

2차 미분은

\begin{aligned}
H_1''&=\frac{-6+12r}{L^2},&
H_2''&=\frac{-4+6r}{L},\\
H_3''&=\frac{6-12r}{L^2},&
H_4''&=\frac{-2+6r}{L}.
\end{aligned}

8.3 Kronecker 및 기울기 조건

\begin{array}{c|cccc|cccc}
&H_1&H_2&H_3&H_4&H_1'&H_2'&H_3'&H_4'\\ \hline
x=0&1&0&0&0&0&1&0&0\\
x=L&0&0&1&0&0&0&0&1
\end{array}

따라서

v(0)=v_1,\quad v'(0)=\theta_{z1},\quad
v(L)=v_2,\quad v'(L)=\theta_{z2},

w(0)=w_1,\quad w'(0)=-\theta_{y1},\quad
w(L)=w_2,\quad w'(L)=-\theta_{y2}.

또한 $H_1+H_3=1$이다. Hermite 보간의 회전 함수 $H_2,H_4$는 길이 차원을 가지므로 네 함수를 단순히 모두 더하는 partition-of-unity 검사는 의미가 없다. 대신 상수 병진과 선형 rigid rotation field의 정확 재현을 검사한다.

8.4 6성분 운동학 보간행렬

중심선 운동학 벡터를

\mathbf q(x)=
\begin{bmatrix}u&v&w&\theta_x&\theta_y&\theta_z\end{bmatrix}^T
=\mathbf N_q(x)\mathbf d_e^l

로 두면

\mathbf N_q=
\begin{bmatrix}
N_1^L&0&0&0&0&0&N_2^L&0&0&0&0&0\\
0&H_1&0&0&0&H_2&0&H_3&0&0&0&H_4\\
0&0&H_1&0&-H_2&0&0&0&H_3&0&-H_4&0\\
0&0&0&N_1^L&0&0&0&0&0&N_2^L&0&0\\
0&0&-H_1'&0&H_2'&0&0&0&-H_3'&0&H_4'&0\\
0&H_1'&0&0&0&H_2'&0&H_3'&0&0&0&H_4'
\end{bmatrix}_{6\times12}.

성분형은

q_i(x)=\sum_{a=1}^{12}N_{q,ia}(x)d_a^l.

9. Kinematics

9.1 3차원 단면 운동학

도심에서 국부 단면좌표가 $(y,z)$인 점의 강체 병진·회전 부분은

\begin{bmatrix}U_x\\U_y\\U_z\end{bmatrix}
=
\begin{bmatrix}u\\v\\w\end{bmatrix}
+
\begin{bmatrix}\theta_x\\\theta_y\\\theta_z\end{bmatrix}
\times
\begin{bmatrix}0\\y\\z\end{bmatrix}
=
\begin{bmatrix}
u+z\theta_y-y\theta_z\\
v-z\theta_x\\
w+y\theta_x
\end{bmatrix}.

Euler–Bernoulli 조건을 대입한 축방향 변위와 변형률은

U_x=u-zw'-yv',

\boxed{
\varepsilon_{xx}(x,y,z)
=\frac{\partial U_x}{\partial x}
=\varepsilon_0+z\kappa_y-y\kappa_z
}

이다.

위 식은 이 요소가 보유하는 6개 단면운동 자유도만 나타내며 비원형 단면의 Saint-Venant warping 변위장은 포함하지 않는다. 축 및 굽힘 변형률은 이 운동학에서 유도하고, 비틀림은 별도의 1차원 에너지/단면력 관계 $T=GJ\kappa_x$로 표현한다.

9.2 일반화 변형률

\boxed{
\boldsymbol\varepsilon=
\begin{bmatrix}
\varepsilon_0\\\kappa_x\\\kappa_y\\\kappa_z
\end{bmatrix}
=
\begin{bmatrix}
u'\\\theta_x'\\\theta_y'\\\theta_z'
\end{bmatrix}
=
\begin{bmatrix}
u'\\\theta_x'\\-w''\\v''
\end{bmatrix}}

9.3 Strain–displacement matrix

\boldsymbol\varepsilon(x)=\mathbf B(x)\mathbf d_e^l

이며 $r=x/L$에서 성분행렬은

\boxed{
\mathbf B=
\begin{bmatrix}
-\frac1L&0&0&0&0&0&\frac1L&0&0&0&0&0\\
0&0&0&-\frac1L&0&0&0&0&0&\frac1L&0&0\\
0&0&\frac{6-12r}{L^2}&0&\frac{-4+6r}{L}&0&
0&0&\frac{-6+12r}{L^2}&0&\frac{-2+6r}{L}&0\\
0&\frac{-6+12r}{L^2}&0&0&0&\frac{-4+6r}{L}&
0&\frac{6-12r}{L^2}&0&0&0&\frac{-2+6r}{L}
\end{bmatrix}_{4\times12}}

즉 비영 성분은

\begin{aligned}
\varepsilon_0&=\frac{u_2-u_1}{L},\\
\kappa_x&=\frac{\theta_{x2}-\theta_{x1}}{L},\\
\kappa_y&=\frac{6-12r}{L^2}w_1
+\frac{-4+6r}{L}\theta_{y1}
+\frac{-6+12r}{L^2}w_2
+\frac{-2+6r}{L}\theta_{y2},\\
\kappa_z&=\frac{-6+12r}{L^2}v_1
+\frac{-4+6r}{L}\theta_{z1}
+\frac{6-12r}{L^2}v_2
+\frac{-2+6r}{L}\theta_{z2}.
\end{aligned}

10. Constitutive Contract

승인된 *BEAM GENERAL SECTION, SECTION=GENERAL 수학 mapping은 입력 순서 A, I11, I12, I22, J에 대해

\boxed{I_y=I11,\qquad I_z=I22,\qquad I12=0}

이다. 여기서 section axis 1은 n1 -> local y, section axis 2는 t x n1 -> local z이므로 $I11$은 local $y$축 굽힘의 I_y, $I22$는 local $z$축 굽힘의 $I_z$에 대응한다. $I12\ne0$을 대각 행렬로 근사하지 않는다.

10.1 단면 적분

도심축 조건은

\int_A y\,dA=0,\qquad \int_A z\,dA=0

이고 주축 조건은

I_{yz}=\int_A yz\,dA=0.

단면 특성은

A=\int_A dA,\qquad
I_y=\int_A z^2dA,\qquad
I_z=\int_A y^2dA.

축응력과 단면력 정의는

\sigma_{xx}=E\varepsilon_{xx},

N=\int_A\sigma_{xx}dA,\qquad
M_y=\int_Az\sigma_{xx}dA,\qquad
M_z=-\int_Ay\sigma_{xx}dA.

이를 적분하면

N=EA\varepsilon_0,\qquad
M_y=EI_y\kappa_y,\qquad
M_z=EI_z\kappa_z.

Saint-Venant 비틀림은

T=GJ\kappa_x

로 둔다. 여기서 $J$는 polar moment와 항상 같은 값이 아니라, 일반 단면에 대한 Saint-Venant torsional constant이다.

10.2 일반화 탄성행렬

\boxed{
\mathbf s=\mathbf D\boldsymbol\varepsilon,
\qquad
\mathbf D=
\begin{bmatrix}
EA&0&0&0\\
0&GJ&0&0\\
0&0&EI_y&0\\
0&0&0&EI_z
\end{bmatrix}_{4\times4}}

성분형은

s_i=\sum_{j=1}^{4}D_{ij}\varepsilon_j,

N=EA\varepsilon_0,\quad
T=GJ\kappa_x,\quad
M_y=EI_y\kappa_y,\quad
M_z=EI_z\kappa_z.

필수 유효성 조건은

E>0,\quad G>0,\quad A>0,\quad J>0,\quad I_y>0,\quad I_z>0.

V0의 균질 등방성 재료는 입력의 정확한 $E,\nu$로부터

G=\frac{E}{2(1+\nu)}

를 계산하고 위의 E>0, G>0 조건을 적용한다. 독립 G 입력은 승인된 V0 입력 계약이 아니다.

10.3 비주축 단면에 대한 참고식

현재 범위 밖이지만 $I_{yz}\ne0$이면 굽힘 블록은

\begin{bmatrix}M_y\\M_z\end{bmatrix}
=E
\begin{bmatrix}
I_y&-I_{yz}\\
-I_{yz}&I_z
\end{bmatrix}
\begin{bmatrix}\kappa_y\\\kappa_z\end{bmatrix}

가 된다. 현재 구현 계약에서 이 항을 0으로 조용히 버리면 안 되며 입력을 거부하거나 향후 별도 기능으로 확장해야 한다.

11. Element Equations

11.1 내부력, 잔차 및 강성

\mathbf f_{\mathrm{int},e}^l
=\int_0^L\mathbf B^T\mathbf s\,dx
=\int_0^L\mathbf B^T\mathbf D\mathbf B\,dx\,\mathbf d_e^l
=\mathbf K_e^l\mathbf d_e^l.

성분형은

f_{\mathrm{int},a}^l
=\int_0^L\sum_{i=1}^{4}B_{ia}s_i\,dx,

K_{ab}^l
=\int_0^L\sum_{i=1}^{4}\sum_{j=1}^{4}
B_{ia}D_{ij}B_{jb}\,dx.

요소 잔차는

\boxed{\mathbf r_e^l=\mathbf K_e^l\mathbf d_e^l-\mathbf f_{\mathrm{ext},e}^l}

이다. 선형 문제이므로 tangent와 stiffness가 같고 변위와 무관하다.

11.2 부분 강성행렬

축 자유도 [u_1,u_2]:

\mathbf K_a=\frac{EA}{L}
\begin{bmatrix}1&-1\\-1&1\end{bmatrix}.

비틀림 자유도 [\theta_{x1},\theta_{x2}]:

\mathbf K_t=\frac{GJ}{L}
\begin{bmatrix}1&-1\\-1&1\end{bmatrix}.

$v$–\theta_z 자유도 [v_1,\theta_{z1},v_2,\theta_{z2}]:

\mathbf K_{bz}=\frac{EI_z}{L^3}
\begin{bmatrix}
12&6L&-12&6L\\
6L&4L^2&-6L&2L^2\\
-12&-6L&12&-6L\\
6L&2L^2&-6L&4L^2
\end{bmatrix}.

$w$–\theta_y 자유도 [w_1,\theta_{y1},w_2,\theta_{y2}]:

\mathbf K_{by}=\frac{EI_y}{L^3}
\begin{bmatrix}
12&-6L&-12&-6L\\
-6L&4L^2&6L&2L^2\\
-12&6L&12&6L\\
-6L&2L^2&6L&4L^2
\end{bmatrix}.

대표 적분 성분은

K_{v_1v_1}=EI_z\int_0^L(H_1'')^2dx=\frac{12EI_z}{L^3},

K_{v_1\theta_{z1}}=EI_z\int_0^LH_1''H_2''dx
=\frac{6EI_z}{L^2},

K_{w_1\theta_{y1}}=EI_y\int_0^L(-H_1'')H_2''dx
=-\frac{6EI_y}{L^2}.

마지막 식의 음수가 $\theta_y=-w'$에서 직접 발생한다.

11.3 전체 12×12 국부 강성행렬

다음 약어를 둔다.

\begin{gathered}
a=\frac{EA}{L},\qquad t=\frac{GJ}{L},\\
b_y=\frac{12EI_y}{L^3},\quad c_y=\frac{6EI_y}{L^2},
\quad d_y=\frac{4EI_y}{L},\quad e_y=\frac{2EI_y}{L},\\
b_z=\frac{12EI_z}{L^3},\quad c_z=\frac{6EI_z}{L^2},
\quad d_z=\frac{4EI_z}{L},\quad e_z=\frac{2EI_z}{L}.
\end{gathered}

그러면 4.2절의 자유도 순서에 대한 국부 강성행렬은

\boxed{
\mathbf K_e^l=
\begin{bmatrix}
a&0&0&0&0&0&-a&0&0&0&0&0\\
0&b_z&0&0&0&c_z&0&-b_z&0&0&0&c_z\\
0&0&b_y&0&-c_y&0&0&0&-b_y&0&-c_y&0\\
0&0&0&t&0&0&0&0&0&-t&0&0\\
0&0&-c_y&0&d_y&0&0&0&c_y&0&e_y&0\\
0&c_z&0&0&0&d_z&0&-c_z&0&0&0&e_z\\
-a&0&0&0&0&0&a&0&0&0&0&0\\
0&-b_z&0&0&0&-c_z&0&b_z&0&0&0&-c_z\\
0&0&-b_y&0&c_y&0&0&0&b_y&0&c_y&0\\
0&0&0&-t&0&0&0&0&0&t&0&0\\
0&0&-c_y&0&e_y&0&0&0&c_y&0&d_y&0\\
0&c_z&0&0&0&e_z&0&-c_z&0&0&0&d_z
\end{bmatrix}}

이어야 한다.

11.4 국부 요소 끝력의 성분형

분포하중이 없는 경우

\mathbf f_{\mathrm{int},e}^l=
\begin{bmatrix}
F_{x1}&F_{y1}&F_{z1}&M_{x1}&M_{y1}&M_{z1}&
F_{x2}&F_{y2}&F_{z2}&M_{x2}&M_{y2}&M_{z2}
\end{bmatrix}^T
=\mathbf K_e^l\mathbf d_e^l.

각 성분은

\begin{aligned}
F_{x1}&=a(u_1-u_2),\\
F_{y1}&=b_z(v_1-v_2)+c_z(\theta_{z1}+\theta_{z2}),\\
F_{z1}&=b_y(w_1-w_2)-c_y(\theta_{y1}+\theta_{y2}),\\
M_{x1}&=t(\theta_{x1}-\theta_{x2}),\\
M_{y1}&=-c_yw_1+d_y\theta_{y1}+c_yw_2+e_y\theta_{y2},\\
M_{z1}&=c_zv_1+d_z\theta_{z1}-c_zv_2+e_z\theta_{z2},\\[2mm]
F_{x2}&=-a(u_1-u_2),\\
F_{y2}&=-b_z(v_1-v_2)-c_z(\theta_{z1}+\theta_{z2}),\\
F_{z2}&=-b_y(w_1-w_2)+c_y(\theta_{y1}+\theta_{y2}),\\
M_{x2}&=-t(\theta_{x1}-\theta_{x2}),\\
M_{y2}&=-c_yw_1+e_y\theta_{y1}+c_yw_2+d_y\theta_{y2},\\
M_{z2}&=c_zv_1+e_z\theta_{z1}-c_zv_2+d_z\theta_{z2}.
\end{aligned}

11.5 차원과 기본 행렬 특성

\mathbf B:4\times12,\quad
\mathbf D:4\times4,\quad
\mathbf K_e^l:12\times12,\quad
\mathbf d_e^l,\mathbf f_e^l:12\times1.

정확한 산술에서

(\mathbf K_e^l)^T=\mathbf K_e^l,
\qquad
(\mathbf d_e^l)^T\mathbf K_e^l\mathbf d_e^l\ge0,
\qquad
\operatorname{rank}(\mathbf K_e^l)=6.

나머지 6개 null modes는 자유 요소의 강체운동이다. 구조 경계조건을 적용한 뒤에만 전역 자유도 강성행렬이 positive definite일 수 있다.

Mass와 damping은 현재 선형 정적 범위에서 N/A이다.

12. External Force Vector

12.1 일반 일관 분포하중

이 절과 12.2절은 요소 formulation을 검증하는 constant local line-load kernel의 수학 계약이다. V0 parser와 CLI는 Abaqus *DLOAD를 지원하거나 Domain distributed-load object를 만들지 않으므로 이 kernel을 해석 pipeline의 지원 하중으로 해석하지 않는다.

국부 단위길이당 하중 벡터를

\mathbf p_l(x)=
\begin{bmatrix}p_x&p_y&p_z&m_x&0&0\end{bmatrix}^T

로 두면 일관 요소하중은

\boxed{
\mathbf f_{\mathrm{dist},e}^l
=\int_0^L\mathbf N_q^T(x)\mathbf p_l(x)\,dx}

이다. 성분형은

f_a^l=\int_0^L\sum_{i=1}^{6}N_{q,ia}(x)p_{l,i}(x)\,dx.

분포 $m_y,m_z$까지 지원하려면 $\mathbf p_l$의 마지막 두 성분과 이에 대응하는 strong form을 별도 I/O 계약에서 명시해야 한다. 현재 baseline에서는 제외한다.

12.2 상수 분포하중의 닫힌형

$p_x,p_y,p_z,m_x$가 요소에서 상수이면

\boxed{
\mathbf f_{\mathrm{dist},e}^l=
\begin{bmatrix}
\frac{p_xL}{2}\\
\frac{p_yL}{2}\\
\frac{p_zL}{2}\\
\frac{m_xL}{2}\\
-\frac{p_zL^2}{12}\\
\frac{p_yL^2}{12}\\
\frac{p_xL}{2}\\
\frac{p_yL}{2}\\
\frac{p_zL}{2}\\
\frac{m_xL}{2}\\
\frac{p_zL^2}{12}\\
-\frac{p_yL^2}{12}
\end{bmatrix}}

이다. 벡터의 5, 11번째 M_y 성분 부호는 w 보간에 $-H_2,-H_4$가 사용되기 때문에 발생한다.

12.3 요소 내부 집중하중과 끝하중

$x_p$에 작용하는 국부 집중 일반화 하중 $\mathbf P_p$의 등가절점하중은

\mathbf f_p^l=\mathbf N_q^T(x_p)\mathbf P_p.

절점 끝하중은 해당 요소 자유도 성분에 직접 더한다. 동일한 물리 하중을 구조 절점 하중과 요소 끝하중으로 중복 조립하지 않아야 한다.

13. Numerical Integration

13.1 강성 적분

자연좌표에서

\mathbf K_e^l
=\int_{-1}^{1}\mathbf B^T(\xi)\mathbf D\mathbf B(\xi)J_x\,d\xi.

상수 물성에 대해 $\mathbf B$의 굽힘행은 $\xi$의 1차식이므로 $\mathbf B^T\mathbf D\mathbf B$는 최대 2차식이다. 2점 Gauss–Legendre 적분은 3차 다항식까지 정확하므로 전체 강성을 정확히 적분한다.

\xi_1=-\frac1{\sqrt3},\qquad
\xi_2=+\frac1{\sqrt3},\qquad
w_1=w_2=1.

행렬 및 성분형은

\boxed{
\mathbf K_e^l
=\sum_{g=1}^{2}\mathbf B^T(\xi_g)\mathbf D\mathbf B(\xi_g)J_xw_g}

K_{ab}^l
=\sum_{g=1}^{2}\sum_{i=1}^{4}\sum_{j=1}^{4}
B_{ia}(\xi_g)D_{ij}B_{jb}(\xi_g)J_xw_g.

1점 적분은 굽힘 부분에 spurious zero-energy mode를 만들 수 있으므로 사용하지 않는다. 구현 초기에는 2점 Gauss 결과를 11.3절 닫힌형 행렬과 성분별 비교한다.

13.2 하중 적분

상수 분포하중에 대한 cubic Hermite 함수 적분은 2점 Gauss로 정확하다.

\mathbf f_{\mathrm{dist},e}^l
=\sum_{g=1}^{2}\mathbf N_q^T(\xi_g)
\mathbf p_l(\xi_g)J_xw_g.

요소 내에서 물성이나 하중이 변하는 기능을 추가할 때는 다항 차수와 discontinuity를 기준으로 적분차수를 다시 정해야 하며, 현재 2점을 일반 규칙으로 오해하면 안 된다.

14. Local-to-Global Transformation

14.1 3×3 회전행렬

전역 벡터를 국부 성분으로 변환하는 행렬을

\boxed{
\mathbf R=
\begin{bmatrix}
\mathbf e_x^T\\
\mathbf e_y^T\\
\mathbf e_z^T
\end{bmatrix}}

로 정의한다.

\mathbf v^l=\mathbf R\mathbf v^g,
\qquad
v_i^l=\sum_{I=1}^{3}R_{iI}v_I^g.

검사 조건은

\mathbf R\mathbf R^T=\mathbf I_3,\qquad
\det\mathbf R=+1.

14.2 12×12 요소 변환행렬

병진과 small-rotation vector는 같은 회전행렬로 변환한다.

\boxed{
\mathbf T=\operatorname{diag}(\mathbf R,\mathbf R,\mathbf R,\mathbf R)}

첫 두 block은 node 1의 병진/회전, 다음 두 block은 node 2의 병진/회전이다.

\mathbf d_e^l=\mathbf T\mathbf d_e^g,
\qquad
d_a^l=\sum_{A=1}^{12}T_{aA}d_A^g.

14.3 전역 요소식

가상일 불변성에서

(\delta\mathbf d_e^g)^T\mathbf f_e^g
=(\delta\mathbf d_e^l)^T\mathbf f_e^l

이므로

\boxed{
\mathbf K_e^g=\mathbf T^T\mathbf K_e^l\mathbf T,
\qquad
\mathbf f_e^g=\mathbf T^T\mathbf f_e^l.}

성분형은

K_{AB}^g
=\sum_{a=1}^{12}\sum_{b=1}^{12}
T_{aA}K_{ab}^lT_{bB},
\qquad
f_A^g=\sum_{a=1}^{12}T_{aA}f_a^l.

15. Global Assembly and Linear Static Solution

요소 전역 자유도를 전체 자유도에서 추출하는 Boolean/scatter 행렬을 $\mathbf A_e$라 하면

\mathbf d_e^g=\mathbf A_e\mathbf d.

조립식은

\boxed{
\mathbf K=\sum_e\mathbf A_e^T\mathbf K_e^g\mathbf A_e,
\qquad
\mathbf F=\sum_e\mathbf A_e^T\mathbf f_e^g+\mathbf F_{\mathrm{nodal}}.}

성분형은

K_{IJ}=\sum_e\sum_{A=1}^{12}\sum_{B=1}^{12}
A_{e,AI}K_{e,AB}^gA_{e,BJ},

F_I=\sum_e\sum_{A=1}^{12}A_{e,AI}f_{e,A}^g+F_{\mathrm{nodal},I}.

자유/구속 자유도로 partition하면

\begin{bmatrix}
\mathbf K_{ff}&\mathbf K_{fc}\\
\mathbf K_{cf}&\mathbf K_{cc}
\end{bmatrix}
\begin{bmatrix}\mathbf d_f\\\mathbf d_c\end{bmatrix}
=
\begin{bmatrix}\mathbf F_f\\\mathbf F_c\end{bmatrix}.

주어진 $\mathbf d_c$에 대해

\boxed{
\mathbf K_{ff}\mathbf d_f
=\mathbf F_f-\mathbf K_{fc}\mathbf d_c}

를 풀고, 구속 자유도 반력은 원래 조립식에서

\boxed{
\mathbf R_c
=\mathbf K_{cf}\mathbf d_f+\mathbf K_{cc}\mathbf d_c-\mathbf F_c}

로 복구한다. 성분형은

R_I=\sum_JK_{IJ}d_J-F_I,\qquad I\in\mathcal C.

16. Output Recovery

16.1 변위

전역 요소변위에서

\mathbf d_e^l=\mathbf T\mathbf d_e^g,
\qquad
\mathbf q(\xi)=\mathbf N_q(\xi)\mathbf d_e^l

로 국부 nodal/section 변위를 복구한다.

16.2 일반화 변형률과 단면력

두 Gauss point $\xi_g=\pm1/\sqrt3$에서 일반화 변형률과 generalized section resultant를 다음과 같이 복구한다.

\boxed{
\boldsymbol\varepsilon(\xi_g)=\mathbf B(\xi_g)\mathbf d_e^l,
\qquad
\mathbf s(\xi_g)=\mathbf D\mathbf B(\xi_g)\mathbf d_e^l.}

성분형은

\varepsilon_i(\xi_g)=\sum_aB_{ia}(\xi_g)d_a^l,
\qquad
s_i(\xi_g)=\sum_jD_{ij}\varepsilon_j(\xi_g).

endpoint section resultant는 별도로 $\xi=-1,+1$에서 $\mathbf D\mathbf B\mathbf d_e^l$을 직접 평가한 local $[N,T,M_y,M_z]$이다. 이는 양의 local x 면을 기준으로 한 section-cut 부호를 가지며 endpoint 외향 법선 부호를 곱하지 않는다. 상수 단면 요소에서 $N,T$는 요소 내 상수이고 $M_y,M_z$는 선형이다. 따라서 Gauss generalized values와 endpoint section resultant는 위치와 row identity가 서로 다른 결과다.

16.3 축응력

단면점 $(y,z)$의 axial S11은

\boxed{
S11(\xi,y,z)=\sigma_{xx}(\xi,y,z)
=E\left[\varepsilon_0(\xi)
+z\kappa_y(\xi)-y\kappa_z(\xi)\right].}

S11은 두 Gauss point와 입력 *SECTION POINTS의 각 단면점에서 평가한다. 입력 section point가 없으면 centroid $(y,z)=(0,0)$을 source=fesa-default로 사용한다. Abaqus section-point 좌표와의 승인 mapping은 x1=y, $x2=z$이다.

Euler–Bernoulli 가정은 횡전단변형을 제거하므로 이 정식화만으로 transverse shear stress 분포를 constitutive recovery하지 않는다. 또한 일반 단면의 비틀림 전단응력 분포는 J 하나만으로 결정되지 않으므로 baseline output은 $T$까지로 제한한다. 원형 단면처럼 별도 승인된 section recovery가 있을 때만 \tau=T\rho/J 형태를 사용할 수 있으나 이는 V0 범위 밖이다. S11은 analytical/unit test와 HDF5 contract로 검증하고 Abaqus beam-stress reference comparison은 명시적 N/A이다.

16.4 요소 끝력

분포하중을 포함하는 요소의 equilibrium end action은

\boxed{
\mathbf f_{\mathrm{end},e}^l
=\mathbf K_e^l\mathbf d_e^l-\mathbf f_{\mathrm{dist},e}^l}

이고 전역 성분은

\mathbf f_{\mathrm{end},e}^g
=\mathbf T^T\mathbf f_{\mathrm{end},e}^l.

각 endpoint row의 local component 순서는 $[F_X,F_Y,F_Z,M_X,M_Y,M_Z]$이다. V0 parser/CLI는 distributed load를 생성하지 않으므로 승인된 해석 pipeline에서는 $\mathbf f_{\mathrm{dist},e}^l=\mathbf0$이다. 이 경우 \xi_1=-1, $\xi_2=+1$에서 외향 법선 부호를 각각 n_{\mathrm{out},1}=-1, $n_{\mathrm{out},2}=+1$로 두면 5.3절의 자연경계 부호와

F_{X,i}=n_{\mathrm{out},i}N(\xi_i),\qquad
M_{X,i}=n_{\mathrm{out},i}T(\xi_i),\qquad
M_{Y,i}=n_{\mathrm{out},i}M_y(\xi_i),\qquad
M_{Z,i}=n_{\mathrm{out},i}M_z(\xi_i),

F_{Y,i}=-n_{\mathrm{out},i}
\left.\frac{dM_z}{dx}\right|_{\xi_i},\qquad
F_{Z,i}=n_{\mathrm{out},i}
\left.\frac{dM_y}{dx}\right|_{\xi_i}

가 일치해야 한다. 따라서 이 V0 조건의 equilibrium end action은 요소 끝에 작용하는 outward action이고, endpoint section resultant는 동일 section을 양의 local x 면으로 자른 값이다. Abaqus node-station force row에는 equilibrium end action이 아니라 section resultant를 동일한 section-cut 부호로 정규화해 비교하며, interior node의 두 endpoint가 승인 tolerance 안에서 먼저 일치해야 한다. 승인된 component mapping은 SF1 -> N, SM1 -> My, SM2 -> Mz, SM3 -> T이고, CSV에 없는 transverse SF2, SF3는 reference comparison 대상이 아니다.

$\mathbf D\mathbf B\mathbf d_e^l$로 얻는 section resultant와 위 equilibrium end action은 분포하중이 있을 때 서로 다른 recovery 표현이다. cubic displacement만으로 uniform transverse load의 quadratic exact moment diagram을 한 요소에서 표현할 수 없으므로 formulation-only line-load kernel에서는 위 V0 endpoint 관계를 적용하거나 두 결과를 같은 row로 동일시하지 않는다.

16.5 반력

반력은 요소 끝력을 임의로 합하는 별도 경로가 아니라 15절의 전체 조립 잔차 $\mathbf K\mathbf d-\mathbf F$에서 구한다. 이렇게 해야 nodal load, element load 및 비영 prescribed displacement가 모두 반영된다.

16.6 Nodal extrapolation

• N,T: 상수이므로 extrapolation 불필요
• M_y,M_z,\kappa_y,\kappa_z: 선형이므로 두 Gauss point 값에서 끝점으로 선형 extrapolation할 수 있으나, 직접 $\mathbf B(\xi=\pm1)$를 평가하는 방법을 기준으로 한다.
• 여러 요소의 nodal stress averaging은 재료/단면 불연속을 숨길 수 있으므로 현재 baseline에서 수행하지 않는다.

17. Algorithm Pseudocode

for each element e:
    read X1, X2, Abaqus first section axis n1, E, nu, A, I11, I12, I22, J
    require I12 == 0
    map a = n1, Iy = I11, Iz = I22, G = E / (2 * (1 + nu))

    DeltaX = X2 - X1
    L = norm(DeltaX)
    coordinate_scale = max(1, norm(X1), norm(X2))
    require L > 1e-12 * coordinate_scale

    ex = DeltaX / L
    ey_trial = a - dot(a, ex) * ex
    require norm(ey_trial) > 1e-12 * max(1, norm(a))
    ey = ey_trial / norm(ey_trial)
    ez = cross(ex, ey)

    R = rows(ex, ey, ez)
    require R * transpose(R) approximately equals identity
    require determinant(R) approximately equals +1
    T = block_diag(R, R, R, R)

    validate E, G, A, Iy, Iz, J are positive
    D = diag(E*A, G*J, E*Iy, E*Iz)

    initialize Kl[12,12] = 0
    initialize fl[12] = 0

    for (xi, weight) in {(-1/sqrt(3), 1), (+1/sqrt(3), 1)}:
        Jx = L / 2
        r = (1 + xi) / 2

        evaluate N1L, N2L
        evaluate H1...H4 and first/second x-derivatives
        form Nq[6,12]
        form B[4,12]

        Kl += transpose(B) * D * B * Jx * weight
        if verifying the formulation-only constant local line-load kernel:
            fl += transpose(Nq) * local_line_load(xi) * Jx * weight

    in verification builds/checks:
        compare Kl component-wise with the closed-form matrix in Section 11.3
        check symmetry, rank 6, and rigid-body null modes

    Kg = transpose(T) * Kl * T
    fg = transpose(T) * fl
    assemble Kg and fg with the element-to-global DOF map

partition global free/constrained DOFs
solve Kff * df = Ff - Kfc * dc
recover Rc = Kcf * df + Kcc * dc - Fc

for each element output recovery:
    de_global = gather global nodal DOFs
    de_local = T * de_global

    for each mandatory Gauss point xi_g:
        epsilon = B(xi_g) * de_local
        section_resultant = D * epsilon
        recover S11 at each input section point or the fesa-default centroid

    for each endpoint xi in {-1, +1}:
        recover endpoint section_resultant = D * B(xi) * de_local
    recover local equilibrium end actions with the endpoint outward-action convention
recover reactions from the assembled residual, not an end-action sum

이 의사코드는 수학적 계산 순서만 정의하며 C++ class, signature 또는 storage API를 정의하지 않는다.

18. Verification and Benchmarks

Finite Element Program Implementation은 element-local 계산, 전역 조립, 해, 결과 복구를 분리하고 새 요소를 patch/benchmark test로 검증할 것을 요구한다. Finite Element Modeling and Convergence Checks는 compatibility, equilibrium, convergence 및 결과 해석을 함께 검사해야 한다는 근거를 제공한다.

18.1 형상함수 단위검사

1. 8.3절의 값/기울기 Kronecker 조건을 machine precision 범위에서 검사한다.
2. N_1^L+N_2^L=1, $H_1+H_3=1$을 검사한다.
3. 상수 $u,\theta_x$를 정확히 재현한다.
4. 선형 $v=c_0+c_1x$에 대해 \theta_z=c_1, $\kappa_z=0$인지 검사한다.
5. 선형 $w=c_0+c_1x$에 대해 \theta_y=-c_1, $\kappa_y=0$인지 검사한다.

18.2 행렬 불변량

승인된 normalized matrix 기준은

\operatorname{normalized\_error}
(\mathbf K_e^l,(\mathbf K_e^l)^T)\le10^{-12},

\operatorname{normalized\_error}
(\mathbf K_{\mathrm{Gauss}}^l,\mathbf K_{\mathrm{closed}}^l)\le10^{-12}

를 검사한다. 자유 요소의 고유치 또는 rank 검사에서는 6개 rigid-body mode와 6개 positive deformation mode가 있어야 한다.

대표 국부 강체모드는

\mathbf r_{T_x}=
\begin{bmatrix}1&0&0&0&0&0&1&0&0&0&0&0\end{bmatrix}^T,

\mathbf r_{T_y}=
\begin{bmatrix}0&1&0&0&0&0&0&1&0&0&0&0\end{bmatrix}^T,

\mathbf r_{T_z}=
\begin{bmatrix}0&0&1&0&0&0&0&0&1&0&0&0\end{bmatrix}^T,

\mathbf r_{R_x}=
\begin{bmatrix}0&0&0&1&0&0&0&0&0&1&0&0\end{bmatrix}^T,

\mathbf r_{R_y}=
\begin{bmatrix}0&0&0&0&1&0&0&0&-L&0&1&0\end{bmatrix}^T,

\mathbf r_{R_z}=
\begin{bmatrix}0&0&0&0&0&1&0&L&0&0&0&1\end{bmatrix}^T.

각 $k$에 대해

\operatorname{normalized\_residual}
(\mathbf K_e^l\mathbf r_k)\le10^{-10}

이어야 한다. Analytical axial, torsion 및 두 굽힘 평면 응답은 relative $10^{-9}$를 적용한다. Norm의 무차원화와 rank scale은 Numerical Review에서 이 승인 기준의 의미를 바꾸지 않는 계산 규약으로 명시한다.

18.3 기본 해석해 benchmark

길이 $L$의 cantilever 하나를 사용하여 다음을 검사한다.

하중	비교량	해석해
끝 축력 P_x	u(L)	P_xL/(EA)
끝 torque M_x	\theta_x(L)	M_xL/(GJ)
끝 횡력 P_y	v(L)	P_yL^3/(3EI_z)
끝 횡력 P_y	\theta_z(L)	P_yL^2/(2EI_z)
끝 횡력 P_z	w(L)	P_zL^3/(3EI_y)
끝 횡력 P_z	\theta_y(L)	-P_zL^2/(2EI_y)
상수 분포하중 p_y	v(L)	p_yL^4/(8EI_z)
상수 분포하중 p_y	\theta_z(L)	p_yL^3/(6EI_z)

한 개의 cubic Hermite 요소는 끝 집중력/모멘트에 대한 cubic 이하 변위장을 정확히 재현한다. 상수 분포 횡하중의 exact displacement는 4차이므로 mesh refinement에 따른 수렴도 함께 확인한다.

18.4 좌표변환 검사

동일한 beam 문제를 전역 $x,y,z$축과 임의의 3차원 방향으로 rigid rotation하여 해석한다. 변위와 힘을 원래 좌표계로 되돌렸을 때 일치해야 한다.

\mathbf T^T\mathbf T=\mathbf I_{12},
\qquad
(\mathbf d_e^g)^T\mathbf K_e^g\mathbf d_e^g
=(\mathbf d_e^l)^T\mathbf K_e^l\mathbf d_e^l.

18.5 평형 검사

\left\|\sum\mathbf F_{\mathrm{applied}}+
\sum\mathbf R\right\|\le\epsilon_F,

\left\|\sum\mathbf M_{\mathrm{applied}}+
\sum\left(\mathbf X\times\mathbf F\right)+
\sum\mathbf M_{\mathrm{reaction}}\right\|\le\epsilon_M.

승인된 B33 reference row는 같은 model, step/frame, quantity, component별로 Abaqus reference 값만 사용해

\operatorname{reference\_scale}=\max_i|r_i|,\qquad
\operatorname{row\_tolerance}=\operatorname{absolute\_floor}
+10^{-6}\operatorname{reference\_scale}

를 적용한다. 승인된 SI bundle의 displacement/rotation floor는 10^{-9}, force/moment floor는 $10^{-3}$이다. Reference 값을 zero-clamp하거나 row를 제거하지 않으며 missing, extra, duplicate, nonfinite, schema/identity mismatch는 수치 판정 전에 실패한다. 이 reference policy는 위 formulation/analytical tolerance를 대체하지 않는다.

19. Numerical Risks

위험	영향	필수 대응/검사
L\approx0	Jacobian singularity, L^{-3} 발산	조립 전 길이 검증
guide vector와 \mathbf e_x 평행	국부 $y,z$축 미정	projection norm 검사 후 입력 오류
E,G,A,I_y,I_z,J\le0	비물리/indefinite stiffness	물성·단면 검증
짧거나 깊은 beam	전단변형 누락으로 과도하게 stiff	Euler–Bernoulli 적용 한계로 기록; 승인되지 않은 자동 warning/error threshold를 두지 않고 B31은 지원 요소로 재해석하지 않음
open section warping	$GJ/L$만으로 비틀림 부족	현재 범위 밖으로 명시
하중이 shear center를 벗어남	굽힘–비틀림 결합 누락	eccentric load를 명시적 force+torque로 환산하거나 후속 기능화
I_y/I_z 또는 stiffness scale 극단값	전역 condition number 악화	단위/스케일 및 solver conditioning 검사
1점 굽힘 적분	spurious mechanism	2점 Gauss 또는 닫힌형 사용
\theta_y=-w' 부호 누락	$w$–M_y 결과 반전	B, K, cantilever sign test
local/global 변환 방향 혼동	회전된 요소 결과 오류	\mathbf d^l=\mathbf T\mathbf d^g 계약 및 energy invariance test
free structure 또는 불충분 구속	singular global matrix	6 rigid modes와 구속조건 진단
nodal averaging	단면 불연속 stress 오염	baseline에서 averaging 금지

Euler–Bernoulli 요소에는 shear strain 자유도가 없으므로 shear locking은 N/A이다. Volumetric locking과 hourglass control도 N/A이지만, 1점 굽힘 적분에 의한 rank deficiency는 별도의 under-integration 위험이다.

20. Evidence and Wiki Gap

20.1 사용한 FEM wiki 근거

• Beam and Frame Finite Elements: Euler–Bernoulli cubic displacement, nodal displacement/rotation, axial·bending·torsion 결합 및 공간좌표 변환
• Isoparametric Finite Elements: 자연좌표, 형상함수, Jacobian, derivative mapping, quadrature 및 distortion/under-integration 위험
• Finite Element Method: strong/variational statement에서 요소방정식, 조립, 해로 이어지는 displacement FEM 흐름
• Finite Element Program Implementation: element-local 계산, 전역 조립, 결과 복구 및 patch/benchmark verification
• Finite Element Modeling and Convergence Checks: compatibility, equilibrium, mesh convergence 및 결과 해석
• A-First-Course-in-the-Finite-Element-Method: beam, frame, grid, space-frame 및 isoparametric formulation의 교재 근거

20.2 근거의 적용 경계

Tier와 provenance가 정리된 research brief는 docs/research/linear-static-3d-euler-beam-research.md에 있다. Wiki와 원출처는 핵심 beam 이론과 Abaqus component 의미를 제공하고, 이 문서의 DOF 순서, \theta_y=-w' 부호, guide-vector 축, 12×12 행렬, line-load vector, 결과 위치와 tolerance는 승인된 project requirement/design과 결합한 FESA 계약이다. 따라서 exact Abaqus B33 내부 알고리즘 동등성이나 full compatibility를 주장하지 않는다.

Research brief의 PARDISO manual 및 deterministic COO-to-CSR 근거 gap은 solver adapter와 assembly 구현 계획의 handoff이며 이 요소의 수학 정식화를 변경하지 않는다. Endpoint sign-normalization gap은 5.3절과 16.4절의 outward end action/positive-face section-cut 관계로 명시했으며 Numerical Review가 독립 재유도해야 한다.

21. Resolved Decisions, Remaining Limitations, and Downstream Handoff

21.1 Resolved decisions and remaining limitations

1. Requirement와 research source는 Metadata의 실제 Step 0–1 산출물로 해결되었다.
2. Orientation과 section mapping은 n1 -> local y, t x n1 -> local z, I_y=I11, I_z=I22, $I12=0$으로 해결되었다.
3. Geometry validation은 3.1절의 scale-aware length/projection 기준으로 해결되었다.
4. Matrix, residual, analytical 및 B33 reference tolerance는 18절의 승인 기준으로 해결되었다.
5. equilibrium end action, endpoint section resultant, Gauss generalized values 및 section-point S11은 16절에서 위치·component·부호가 구분되었다. Stress reference comparison은 N/A이다.
6. Euler–Bernoulli slenderness는 적용 limitation으로 남는다. 공식 자료의 약 1/15 지침을 FESA warning/error threshold로 바꾸지 않는다.
7. Torsional/transverse shear stress, warping, I12\ne0, B31/Timoshenko 및 *DLOAD parser support는 승인 범위 밖이다.

현재 Numerical Review를 막는 미해결 수학 입력은 없다. ready-for-numerical-review는 독립 수치 검토 준비 상태이며 정식화 최종 승인이나 release readiness를 뜻하지 않는다.

21.2 Numerical Review Agent

• $\theta_y=-w'$와 M_y=-EI_yw'' 부호를 독립적으로 재유도한다.
• 4×12 \mathbf B 적분이 12×12 닫힌형과 일치하는지 symbolic/numerical 검사한다.
• 6 rigid-body modes, rank 6, symmetry 및 energy invariance를 검사한다.
• 자연경계력 $F_y,F_z,M_y,M_z$와 요소 끝력 부호를 검사한다.
• 16.4절의 outward equilibrium end action과 positive-face endpoint section resultant 관계를 양 endpoint에서 독립 재유도한다.
• 2점 Gauss exactness와 1점 적분 rank deficiency를 확인한다.

21.3 I/O Definition Agent

• node coordinates, n1, E,\nu,A,I11,I12,I22,J 입력과 승인 mapping을 정의한다.
• guide vector 평행, zero length, nonpositive property validation을 정의한다.
• local/global displacement, reaction, equilibrium end action, endpoint section resultant, Gauss generalized strain/resultant, section-point axial S11의 HDF5 schema와 위치를 구분한다.
• Abaqus keyword subset은 별도 승인 범위만 정의하며 full compatibility를 주장하지 않는다.

21.4 Reference Model Agent

• axial cantilever, torsion cantilever, y/z 양방향 bending cantilever를 분리한다.
• local 축과 일치하는 모델 및 임의 3차원 회전 모델을 모두 준비한다.
• formulation-only uniform line-load kernel은 unit/analytical coverage로 두고 *DLOAD reference support로 확대하지 않는다.
• 승인된 reference/cantilever beam/의 exact legacy path를 read-only로 유지하고 reference artifact를 생성·수정·복원하지 않는다.

21.5 Implementation Planning Agent

TDD 순서는 다음이 적합하다.

1. Hermite 값/미분 및 mapping 단위검사
2. \mathbf B 성분과 닫힌형 \mathbf K_e^l 비교검사
3. rigid modes, symmetry, rank 및 positive-energy 검사
4. 좌표변환 orthogonality와 rotated-element invariance 검사
5. consistent uniform line-load vector 검사
6. axial/torsion/two-plane cantilever 해석해 검사
7. 조립, prescribed displacement 및 reaction 통합검사
8. 결과 복구와 end-action sign 검사

Production C++ 파일을 변경할 때는 관련 C++ test를 먼저 작성해 RED 실패를 확인한 뒤 GREEN 구현과 MSVC x64 Debug/CTest VERIFY를 수행해야 한다.