FESA/docs/formulation/timoshenko-beam-3d.md

2절점 3D Timoshenko Beam 정식화

1. 범위와 가정

FESA Phase 1 Beam은 소변형, 선형 탄성의 2절점 직선 isoparametric Timoshenko 요소다. 단면은 도심과 전단중심이 일치하는 주축 단면이며 (I_{yz}=0)이다. 단면 offset, warping, 기하·재료 비선형은 포함하지 않는다. FESA는 단위 변환을 하지 않으므로 모든 입력은 하나의 일관 단위계를 사용해야 한다.

정식화의 유한요소 이산화와 수치적분은 Bathe, Finite Element Procedures, 2판과 Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis를 따른다. 3차원 isoparametric Beam의 국부 좌표계는 Bathe와 Bolourchi(1979)의 선형화된 부분을 사용한다. 선택적 감차적분은 Hughes, Taylor와 Kanoknukulchai(1977)의 원칙을 Beam 전단항에 적용한다. Plate 운동학은 사용하지 않는다.

2. 국부 좌표계와 자유도

절점 1에서 절점 2로 향하는 단위벡터를 (\mathbf e_x)로 둔다. 입력 orientation 벡터에서 (\mathbf e_x) 성분을 Gram-Schmidt로 제거하고 정규화한 벡터를 (\mathbf e_y)로 두며,

[ \mathbf e_z=\mathbf e_x\times\mathbf e_y ]

로 정의한다. 따라서 ((\mathbf e_x,\mathbf e_y,\mathbf e_z))는 오른손 직교 기저다. 길이가 0인 요소, 길이가 0인 orientation, 요소축과 평행한 orientation은 오류다.

국부 요소 자유도 벡터의 순서는 다음과 같다.

[ \mathbf q_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}. ]

(u,v,w)는 각각 국부 (x,y,z) 병진이고 (\theta_x,\theta_y,\theta_z)는 오른손 규칙의 회전 성분이다. 전역 자유도도 절점별 ([u_x,u_y,u_z,r_x,r_y,r_z]) 순서를 사용한다.

기준축에서 떨어진 점 ((0,y,z))의 선형화된 변위장은

[ U_x=u+z\theta_y-y\theta_z,\qquad U_y=v-z\theta_x,\qquad U_z=w+y\theta_x ]

이다.

3. 변형률과 부호

일반화 변형률과 대응 단면력의 component 순서는

[ \boldsymbol\varepsilon_s = \begin{bmatrix} \epsilon&\gamma_y&\gamma_z&\kappa_x&\kappa_y&\kappa_z \end{bmatrix}^{T}, \qquad \mathbf s = \begin{bmatrix} N&V_y&V_z&T&M_y&M_z \end{bmatrix}^{T} ]

로 고정한다. 위 변위장과 오른손 회전 부호로부터

[ \begin{aligned} \epsilon &= \frac{du}{dx},& \gamma_y &= \frac{dv}{dx}-\theta_z,& \gamma_z &= \frac{dw}{dx}+\theta_y,\ \kappa_x &= \frac{d\theta_x}{dx},& \kappa_y &= \frac{d\theta_y}{dx},& \kappa_z &= \frac{d\theta_z}{dx} \end{aligned} ]

를 사용한다. 양의 (N,V_y,V_z,T,M_y,M_z)는 위 일반화 변형률과 양의 내부 일률 (\delta\boldsymbol\varepsilon_s^T\mathbf s)을 이루는 방향이다. 이 부호에서 단면 축응력은

[ \sigma_{xx}=E\left(\epsilon+z\kappa_y-y\kappa_z\right) ]

이고 (M_y=\int_A z\sigma_{xx},dA), (M_z=-\int_A y\sigma_{xx},dA)이다.

4. 재료와 단면 constitutive matrix

등방성 선형 탄성의 전단계수는 입력 (E,\nu)로부터

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

로 계산한다. 일반화 constitutive 관계는

[ \mathbf s=\mathbf D\boldsymbol\varepsilon_s,\qquad \mathbf D= \operatorname{diag} \left(EA,;GA_{sy},;GA_{sz},;GJ,;EI_y,;EI_z\right). ]

(A,I_y,I_z,J,A_{sy},A_{sz})는 양의 값이어야 한다. 명시적인 전단강성이 없는 입력은 semantic mapper가 커널 호출 전에 (A_{sy}=A_{sz}=5A/6), SCF=0으로 정규화한다. 커널은 Abaqus의 slenderness compensation을 추가하지 않는다.

5. 자연좌표, 형상함수와 Jacobian

자연좌표는 (\xi\in[-1,1])이고 두 절점의 선형 형상함수는

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

모든 여섯 component를 같은 형상함수로 보간한다. 직선 요소 길이를 (L)이라 하면

[ J=\frac{dx}{d\xi}=\frac{L}{2},\qquad \frac{dN_a}{dx}=\frac{1}{J}\frac{dN_a}{d\xi} ]

이다.

절점 (a)의 자유도 순서 ([u_a,v_a,w_a,\theta_{xa},\theta_{ya},\theta_{za}])에 대한 strain-displacement block은

[ \mathbf B_a(\xi)= \begin{bmatrix} N_{a,x}&0&0&0&0&0\ 0&N_{a,x}&0&0&0&-N_a\ 0&0&N_{a,x}&0&N_a&0\ 0&0&0&N_{a,x}&0&0\ 0&0&0&0&N_{a,x}&0\ 0&0&0&0&0&N_{a,x} \end{bmatrix}, \qquad \mathbf B=\begin{bmatrix}\mathbf B_1&\mathbf B_2\end{bmatrix}. ]

6. 선택적 감차적분과 요소 강성

Constitutive matrix를 축·비틀림·굽힘 부분과 전단 부분으로 나눈다.

[ \begin{aligned} \mathbf D_{ab}&= \operatorname{diag}(EA,0,0,GJ,EI_y,EI_z),\ \mathbf D_s&= \operatorname{diag}(0,GA_{sy},GA_{sz},0,0,0). \end{aligned} ]

국부 강성은

[ \mathbf K_l = \sum_{g=1}^{2} \mathbf B(\xi_g)^T\mathbf D_{ab}\mathbf B(\xi_g)Jw_g + \sum_{g=1}^{1} \mathbf B(\xi_g)^T\mathbf D_s\mathbf B(\xi_g)Jw_g ]

로 직접 적분한다. 축·비틀림·굽힘은 (\xi=\pm1/\sqrt{3},w=1)인 2점 Gauss rule을 사용하고, 전단은 (\xi=0,w=2)인 1점 Gauss rule만 사용한다. 전단항을 2점 적분하지 않는다. 선형 회전장에서 1점 전단 적분은 요소 중앙의 평균 회전을 사용하므로 일정곡률 굽힘 mode의 불필요한 전단에너지를 제거한다.

Matrix12는 matrix[row][column]의 12행 12열 고정 크기 저장을 사용한다. 적분은 각 Gauss point에서 (\mathbf B^T\mathbf D\mathbf B)의 상삼각 entry를 고정된 부동소수점 연산 순서로 누적한 뒤 하삼각에 복사한다.

7. 전역 변환

기존 beam_transformation이 만드는 (\mathbf T)는 각 병진·회전 3-vector를 전역 성분에서 국부 성분으로 투영한다.

[ \mathbf q_l=\mathbf T\mathbf q_g. ]

따라서 전역 강성과 에너지 관계는

[ \mathbf K_g=\mathbf T^T\mathbf K_l\mathbf T,\qquad \frac12\mathbf q_g^T\mathbf K_g\mathbf q_g =\frac12\mathbf q_l^T\mathbf K_l\mathbf q_l ]

이다.

8. 구현과 검증 invariant

커널은 위 적분식으로 local stiffness를 계산하고 변환식으로 global stiffness를 계산한다. 다음 invariant가 정식화 검증 기준이다.

• local/global stiffness는 roundoff 이내에서 대칭이다.
• 세 병진과 세 회전 강체 mode의 변형에너지는 0이다.
• 축 및 비틀림 submatrix는 각각 (EA/L), (GJ/L)을 재현한다.
• 일정곡률 단축·이축 굽힘 에너지는 각각 (EI_y\kappa_y^2L/2), (EI_z\kappa_z^2L/2)이다.
• 일정전단 mode의 에너지는 각각 (GA_{sy}\gamma_y^2L/2), (GA_{sz}\gamma_z^2L/2)이다.
• 일정곡률 mode의 전단에너지는 slenderness 변화와 무관하게 0이며 전체 에너지는 해석값과 일치한다.
• 요소와 자유도를 같은 강체 회전으로 회전하면 변형에너지가 보존된다.

9. 참고문헌과 FESA 차이

• K. J. Bathe, Finite Element Procedures, 2nd ed.
• T. J. R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis.
• K. J. Bathe and S. Bolourchi, “Large Displacement Analysis of Three-Dimensional Beam Structures,” 1979. FESA는 이 문헌의 선형화된 국부 좌표와 변환만 사용한다.
• T. J. R. Hughes, R. L. Taylor, and W. Kanoknukulchai, “A Simple and Efficient Finite Element for Plate Bending,” 1977. FESA는 selective integration 원칙만 사용한다.
• Abaqus 2024, “Choosing a Beam Element” 및 “BEAM GENERAL SECTION.” FESA의 (5A/6) 생략값과 SCF=0은 Phase 1 명시적 계약이며 임의 일반 단면의 보편적인 Abaqus 기본값이 아니다. Abaqus B31 slenderness compensation도 구현하지 않는다.