From 212b6b87476572ae28e642344279011f54855f03 Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Fri, 7 Aug 2026 23:46:44 +0900 Subject: [PATCH] docs: formulate 3D Euler beam element --- ...3d-isoparametric-euler-beam-formulation.md | 1552 +++++++++++++++++ 1 file changed, 1552 insertions(+) create mode 100644 docs/formulations/3d-isoparametric-euler-beam-formulation.md diff --git a/docs/formulations/3d-isoparametric-euler-beam-formulation.md b/docs/formulations/3d-isoparametric-euler-beam-formulation.md new file mode 100644 index 0000000..de791ce --- /dev/null +++ b/docs/formulations/3d-isoparametric-euler-beam-formulation.md @@ -0,0 +1,1552 @@ +# 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|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 + +```text +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|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를 수행해야 한다.