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

3D Isoparametric Euler–Bernoulli Beam Formulation

Metadata

• feature_id: 3d-isoparametric-euler-beam
• source_requirement: 사용자 요청 및 2026-08-07 범위 확인; docs/requirements/3d-isoparametric-euler-beam.md는 아직 없음
• source_research: FEM theory wiki의 Beam and Frame Finite Elements, Isoparametric Finite Elements, A-First-Course-in-the-Finite-Element-Method; docs/research/3d-isoparametric-euler-beam-research.md는 아직 없음
• status: needs-research
• owner_agent: formulation-agent
• date: 2026-08-07

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 요소
• small displacement 및 small rotation
• 선형 정적 해석
• 균질 선형 탄성 재료
• 축 변형, 두 방향 Euler–Bernoulli 굽힘, Saint-Venant 비틀림
• 절점당 자유도 [u,v,w,\theta_x,\theta_y,\theta_z]
• 단면 도심을 지나는 국부 $x$축
• 단면 주축과 국부 $y,z$축의 일치: I_{yz}=0
• 보존력과 일관 등가절점하중
• 자연좌표 $\xi\in[-1,1]$을 사용하는 1차원 등매개 기하 mapping

2.2 제외 범위

• Timoshenko 전단 변형과 shear correction factor
• 기하비선형, 재료비선형, 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 단위

특정 단위계를 고정하지 않고 일관 단위계를 사용한다.

양	차원
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$라 하고 사용자가 제공하는 국부 $y$축 guide vector를 $\mathbf a$라 한다.

\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}.

다음 조건을 만족해야 한다.

L>L_{\mathrm{tol}},\qquad
\|\widetilde{\mathbf e}_y\|>a_{\mathrm{tol}},\qquad
\mathbf e_i\cdot\mathbf e_j=\delta_{ij},\qquad
\mathbf e_x\times\mathbf e_y=\mathbf e_z.

guide vector가 요소 축과 평행하면 국부 $y,z$축과 $I_y,I_z$의 의미가 정해지지 않으므로 입력 오류로 처리한다. 임의의 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}

요소는 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

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.

등방성 재료에서 $G$를 $E,\nu$로 계산한다면

G=\frac{E}{2(1+\nu)},\qquad -1<\nu<\frac12

를 사용한다. $G$가 독립 입력이면 이 변환은 적용하지 않는다.

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 일반 일관 분포하중

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

\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 일반화 변형률과 단면력

요청 위치 $\xi_o$에서

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

성분형은

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

상수 단면 요소에서 $N,T$는 요소 내 상수이고 $M_y,M_z$는 선형이다. 기본 출력 위치는 두 Gauss point이며 필요하면 $\xi=-1,0,+1$에서도 직접 평가한다.

16.3 축응력

단면점 $(y,z)$의 축응력은

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

Euler–Bernoulli 가정은 횡전단변형을 제거하므로 이 정식화만으로 transverse shear stress 분포를 constitutive recovery하지 않는다. 또한 일반 단면의 비틀림 전단응력 분포는 J 하나만으로 결정되지 않으므로 baseline output은 $T$까지로 제한한다. 원형 단면처럼 별도 승인된 section recovery가 있을 때만 \tau=T\rho/J 형태를 사용한다.

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.

$\mathbf D\mathbf B\mathbf d_e^l$로 얻는 section resultant와 위 equilibrium end action은 분포하중이 있을 때 서로 다른 recovery 표현이다. cubic displacement만으로 uniform transverse load의 quadratic exact moment diagram을 한 요소에서 표현할 수 없으므로 두 결과를 같은 데이터로 혼동하지 않는다.

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, local-y guide a, E, G, A, Iy, Iz, J

    DeltaX = X2 - X1
    L = norm(DeltaX)
    require L > length_tolerance

    ex = DeltaX / L
    ey_trial = a - dot(a, ex) * ex
    require norm(ey_trial) > orientation_tolerance
    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
        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 requested element output location xi_o:
    de_global = gather global nodal DOFs
    de_local = T * de_global
    epsilon = B(xi_o) * de_local
    section_resultant = D * epsilon
    recover sigma_xx at requested section points
    recover local/global equilibrium end actions

이 의사코드는 수학적 계산 순서만 정의하며 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 행렬 불변량

허용오차 $\epsilon_K$에 대해

\|\mathbf K_e^l-(\mathbf K_e^l)^T\|\le\epsilon_K,

\|\mathbf K_{\mathrm{Gauss}}^l-\mathbf K_{\mathrm{closed}}^l\|
\le\epsilon_K

를 검사한다. 자유 요소의 고유치 또는 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$에 대해

\|\mathbf K_e^l\mathbf r_k\|\le\epsilon_{\mathrm{rb}}

이어야 한다.

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.

절대/상대/norm tolerance 값은 Requirement 및 Numerical Review 단계에서 별도로 승인해야 한다.

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	slenderness 적용범위 명시; Timoshenko 요소로 전환
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 Wiki gap

현재 wiki는 핵심 이론과 source provenance는 제공하지만, 이 문서에서 고정한 [u,v,w,\theta_x,\theta_y,\theta_z] 순서, \theta_y=-w' 부호, guide-vector 좌표계, 전체 12×12 행렬, 일관 line-load vector 및 HDF5 후보 결과 위치를 한 페이지에 모은 구현 수준 계약은 제공하지 않는다. 이 세부식은 위 근거와 사용자가 승인한 범위를 결합해 유도했으며 Numerical Review를 거쳐야 한다.

21. Open Issues and Downstream Handoff

21.1 공통 open issues

1. 승인된 요구조건 문서 docs/requirements/3d-isoparametric-euler-beam.md가 아직 없다.
2. Tier가 부여된 연구 브리프 docs/research/3d-isoparametric-euler-beam-research.md가 아직 없다.
3. Euler–Bernoulli 적용 slenderness 기준과 경고/오류 정책이 미정이다.
4. absolute, relative 및 norm tolerance가 미정이다.
5. 단면 orientation 입력 규약과 J,I_y,I_z 제공 방식이 I/O 단계에서 미정이다.
6. section resultant와 equilibrium end action의 HDF5 dataset 구분이 미정이다.
7. 일반 단면의 torsional shear stress recovery는 현재 범위 밖이다.

이 항목들 때문에 문서 상태는 needs-research이며, 정식화 자체를 최종 승인한 것으로 간주하지 않는다.

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$와 요소 끝력 부호를 검사한다.
• 2점 Gauss exactness와 1점 적분 rank deficiency를 확인한다.

21.3 I/O Definition Agent

• node coordinates, section local-y guide, E,G,A,I_y,I_z,J 입력을 정의한다.
• guide vector 평행, zero length, nonpositive property validation을 정의한다.
• local/global displacement, reaction, end action, section resultant, generalized strain, axial stress의 HDF5 schema와 위치를 구분한다.
• Abaqus keyword subset은 별도 승인 범위만 정의하며 full compatibility를 주장하지 않는다.

21.4 Reference Model Agent

• axial cantilever, torsion cantilever, y/z 양방향 bending cantilever를 분리한다.
• local 축과 일치하는 모델 및 임의 3차원 회전 모델을 모두 준비한다.
• uniform line load, end force, end moment 및 mixed 3D load 모델을 포함한다.
• reference artifact 생성이나 수정은 승인된 reference phase에서만 수행한다.

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를 수행해야 한다.