FESADev/docs/mitc4-shell-feast-eigen/formulation-mass.md

MITC4 consistent mass와 nodal block lumping 유도

• Feature ID: mitc4-shell-feast-eigen
• Owner: formulation-agent
• Input contract revision: r1
• 기본 보고서: formulation.md

근거는 research evidence의 S01, S02, P01, R01~R03이다. 아래 유도는 FESA의 승인된 운동학에 대한 것이며 Abaqus 내부 질량과의 동등성 주장이 아니다.

M1. 변수, 원형상과 운동학

길이·질량·시간 차원을 각각 L, M, T라 한다. 각 element는 원형상 좌표 X_i [L], 기존 patch에서 얻은 공통 unit nodal director d_i, 일정한 centered thickness t>0 [L], section의 material에서 결정된 일정 density rho>0 [M/L³]를 사용한다. E>0 [M/(LT²)], -1<nu<0.5와 기존 단일층 등방성 plane-stress 재료를 유지한다. 모든 입력 및 중간 결과는 finite여야 한다. 여기의 t는 두께이며 시간은 tau로 쓴다.

자연 좌표는 (xi,eta,zeta) in [-1,1]^3, source node 순서는 (-1,-1),(1,-1),(1,1),(-1,1)이다. 다음 형상함수와 실제 volume mapping을 사용한다.

[ N_i=\tfrac14(1+\xi_i\xi)(1+\eta_i\eta),\quad z=t\zeta/2,\quad X=\sum_iN_i(X_i+z d_i),\quad J=\det[X_{,\xi},X_{,\eta},X_{,\zeta}], \quad dV=J,d\xi d\eta d\zeta. \tag{M1} ]

sum Ni=1, Ni>=0이다. Interpolated director를 normalization하여 M1의 geometry를 변경하지 않는다. 기존 kernel의 EvaluateGeometry와 같은 Jacobian 및 부호를 사용한다. 기존 정식화의 §§35, 811이 director 생성, tangent frame (a_i,b_i,d_i), 기하 허용 조건과 MITC strain을 소유한다.

전역 nodal 좌표는 q_i=[u_X,u_Y,u_Z,theta_X,theta_Y,theta_Z]^T이다. u [L], theta [1]이고 양의 director 변화는 theta cross d이다. 다음 map을 정의한다.

[ S_i=-[d_i]_\times,\quad P_i=I-d_i d_i^T=S_i^TS_i,\quad G_i(z)=[I_3\ \ zS_i],\quad H=[N_1G_1\ \cdots\ N_4G_4]. \tag{M2} ]

따라서 입자 변위/속도는 u(X)=Hq, v(X)=H q_dot이다. Physical nodal coordinates [u_X,u_Y,u_Z,alpha,beta]에서 회전 열은 [-z b_i, z a_i]이다. alpha=a_i^T theta, beta=b_i^T theta, gamma=d_i^T theta이며 S_i d_i=0이다. 20↔24 변환은 기존 PhysicalTransformation20/DrillingTransformation4 그대로이다. MITC assumed shear와 5/6 shear factor는 strain/stiffness에만 작용한다. 속도 보간과 질량에 이 factor를 적용하지 않는다.

M2. Strong/weak form과 기존 stiffness

원형상·무응력 상태의 소변형 무감쇠 자유진동에 대해 continuum parent는 rho u_ddot = div sigma, sigma=C_ps:epsilon, sigma_33=0이다. 자유 표면·자유 edge의 generalized traction은 0이고 essential data는 동차 global DOF 제약이다. 허용한 thickness 운동학과 MITC projection에 제한한 약형은

[ \int_{\Omega_0}\rho,\delta u^T\ddot u,dV+ \int_{\Omega_0}\delta\bar e^T C_5\bar e,dV+ \delta\gamma^T K_d\gamma=0. \tag{M3} ]

원형상 elastic Hessian은 기존 stiffness다.

[ C_5=\operatorname{diag}\left( \frac{E}{1-\nu^2}\begin{bmatrix}1&\nu&0\\nu&1&0\0&0&(1-\nu)/2\end{bmatrix}, \frac56\frac{E}{2(1+\nu)}I_2\right), \quad K_{20}=\sum_{g\in2\times2\times2}\bar B_g^TC_5\bar B_gJ_gw_g, ] [ K_e=T_p^TK_{20}T_p+T_d^T(k_d I_4)T_d,\quad k_d=10^{-3}\min{(K_{20})_{rr}>0:r\text{는 tangent rotation}}. \tag{M4} ]

bar e=[epsilon11,epsilon22,gamma12,gamma13,gamma23]이고 M4는 새 선택이 아니라 기존 정식화 §12 및 실제 kernel의 고정 계약이다. Drilling energy는 수치 stabilization이며 물리 strain/stress가 아니다. K_e q_e는 내부 elastic generalized force, M_L,e q_ddot,e는 inertia다. 최종 반이산식은 M_L q_ddot + Kq=0이다. Harmonic ansatz q(phi,tau)=phi cos(omega tau)에서 Kphi=lambda M_L phi, lambda=omega² [T^-2]를 얻는다. 새 constitutive/history state, geometric stiffness, damping, preload, nonlinear tangent, time integrator 또는 modal stress/reaction recovery는 없다.

M3. Consistent matrix와 충분한 quadrature

T_C=1/2 integral rho |H q_dot|² dV에서

[ C^{(k)}{ij}=\int{\Omega_0}\rho N_iN_j z^k dV\quad(k=0,1,2), \qquad (M_C){ij}=\begin{bmatrix} C^{(0)}{ij}I_3&C^{(1)}{ij}S_j\ C^{(1)}{ij}S_i^T&C^{(2)}_{ij}S_i^TS_j \end{bmatrix}. \tag{M5} ]

세 scalar moment matrix는 각각 symmetric 4×4이다. 차원은 C(k) [M L^k], mass의 uu/u-theta/theta-theta block은 각각 M/ML/ML²이다. M_C [24×24]는 global block 형식이며 20 physical mass의 congruence와 동일하다. 원본 M_C와 그것을 구성한 C(0..2)를 완성·검증한 후 lumping 입력으로 사용한다. Global consistent sparse matrix를 eigen solve에 조립할 필요는 없지만 계산 원본을 lumped 결과로 덮어쓰거나 별도의 mass를 생성한 뒤 consistent에서 유도했다고 부르지 않는다.

M1은 각 자연 좌표에 최대 1차이다. Jacobian의 세 열의 최대 다항 차수는 각각 (0,1,1), (1,0,1), (1,1,0)이므로 determinant의 성분별 최대 차수는 (2,2,2)다. Ni Nj z^k J의 차수는 최대 (4,4,k+2), 즉 (4,4,4)이다. 따라서 질량에는 3×3×3 Gauss–Legendre를 사용한다. 각 축의 points는 [-sqrt(3/5),0,+sqrt(3/5)], weights는 [5/9,8/9,5/9]이다. 3-point rule의 5차 정확성으로 M5의 모든 항을 이 geometry에서 정확 적분한다 (floating-point roundoff 제외). 추가 density/thickness 변화를 가정한 결론이 아니다.

Loop 순서는 증가하는 xi, 그 안의 eta, 가장 안의 zeta, 이후 i<=j, k=0,1,2로 고정하고 대칭 항을 복사한다. Mass 계산의 27개 실제 평가점에서도 기존 geometry predicate를 적용한다. J<=0, nonfinite frame/basis 등의 실패를 weight 절댓값, 점 생략이나 flat replacement로 수리하지 않는다. 기존 stiffness 8-point 경로는 변경하지 않는다. 이 검사는 27점에서의 positivity를 보장하며 volume 모든 위치의 injectivity에 대한 새 theorem이나 기존 warp/curvature 범위의 확대를 주장하지 않는다.

왜곡된 flat element라도 같은 식을 쓴다. 서로 다른 director가 있으면 J(zeta)의 1차 항 때문에 C(1)이 0이 아닐 수 있다. Centered section만을 이유로 coupling을 버리거나 C(2)=t² C(0)/12를 일반식으로 사용하지 않는다.

M4. 계산된 consistent moment에서 nodal block 유도

각 node의 세 row sum과 질량 block을 정의한다.

[ m_i=\sum_j C^{(0)}{ij},\quad h_i=\sum_j C^{(1)}{ij},\quad j_i=\sum_j C^{(2)}{ij},\qquad L_i=\begin{bmatrix}m_iI_3&h_iS_i\h_iS_i^T&j_iP_i\end{bmatrix}, \quad M{L,e}=\operatorname{blockdiag}(L_1,L_2,L_3,L_4). \tag{M6} ]

이는 scalar moment row sum 후 tensor로 재구성하는 nodal block lumping이다. Global 전체 행의 합, absolute row sum, diagonal scaling 또는 nodal quadrature와 혼용하지 않는다. 같은 source node의 incident element block은 기존 source element 순서로 합한다. 공통 d_i 때문에 조립 후에도 M6 형태다. Local와 global에서 diagonal이 아닌 6×6 coupling은 허용한 lumped mass의 일부다.

sum Nj=1을 M5에 대입하면 정확히

[ (m_i,h_i,j_i)=\int\rho N_i(1,z,z^2)dV,\quad L_i=\int\rho N_i G_i(z)^TG_i(z)dV. \tag{M7} ]

m_i>0, j_i>0, m_i j_i-h_i²>0이다. 마지막 부등식은 양의 measure에 대해 z가 상수가 아닌 점(양의 두께 및 서로 다른 thickness Gauss점)과 strict Cauchy–Schwarz에서 따른다. 이를 이용하면

[ q_i^TL_iq_i=m_i\left|u_i+\frac{h_i}{m_i}S_i\theta_i\right|^2+ \left(j_i-\frac{h_i^2}{m_i}\right)|S_i\theta_i|^2. \tag{M8} ]

따라서 L_i의 rank는 5, null은 정확히 [0;d_i]의 span이다. zbar_i=h_i/m_i [L], j_i-h_i²/m_i [ML²]는 node에 할당된 thickness measure의 centroid와 central tangent inertia다. 음/영 variance나 nonfinite를 clamp하지 않고 mass 검증 실패로 처리한다. 모든 노드에 최소 한 활성 shell이 연결되어야 한다.

M6의 유도 연결은 C(k) 원본→row sum→block reconstruction으로 추적한다. 동시에 M7을 독립 직접 적분하여 두 결과의 일치를 검증할 수 있다. 서로 다른 nodal frame의 rotary component를 직접 더하지 않는다.

M5. 보존되는 양과 근사되는 양

모든 등식은 동일 원형상·density와 quadrature measure에 관한 것이다.

1. 총 병진 질량 m=sum mi=integral rho dV 및 모든 uniform translation의 kinetic energy가 consistent와 정확히 같다. 각 방향의 total mass가 m이며 trace(Muu)를 총 질량이라고 정의하지 않는다.
2. Thickness scalar moment는 sum hi=integral rho z dV, sum ji=integral rho z² dV다. 각 nodal tangent tensor는 ji Pi이고 drilling inertia는 0이다. 이는 물리 director의 회전 관성 정의이며 모든 rigid rotation의 전체 tensor를 consistent와 같게 만드는 계약은 아니다.
3. Lumping measure의 위치를 x_i(z)=X_i+z d_i로 보면 spatial first moment sum_i(mi Xi+hi di)=integral rho X dV이므로 center of mass가 보존된다. Node의 Xi에만 mi를 놓고 hi di를 누락하면 곡면에서 이 등식이 깨진다.
4. v_i(z)=u_dot_i+z S_i theta_dot_i에 대해 Jensen identity가 다음을 준다.

[ \dot q^T(M_L-M_C)\dot q= \frac12\int\rho\sum_{i,j}N_iN_j|v_i-v_j|^2dV\ \ge0. \tag{M9} ]

두 kinetic energy의 차이는 M9의 절반이다. Equality는 positive-weight point마다 모든 v_i가 같을 때 성립한다. 일반 속도장 및 rigid angular motion의 orbital inertia는 달라질 수 있으며 lumped 값이 그 quadratic form에서 작아지지 않는다. Rigid spatial rotation은 u_i=Omega cross Xi, theta_i=Pi Omega로 재현된다. 이때 M9는 nodal 위치 분산으로 생기는 angular-inertia 차이를 정확히 설명한다. 두 spectrum의 동일성, arbitrary mesh의 특정 수렴 차수 또는 Abaqus lumping 동등성을 acceptance로 요구하지 않는다. Smooth refinement에서 차이가 감소하는지는 독립 evidence다.

Proper orthogonal global rotation Q에 대해 D_Q=blockdiag(Q,Q)를 두면 S(Qd)=Q S(d) Q^T, L'_i=D_Q L_i D_Q^T이며 consistent에도 같은 congruence가 성립한다. 이는 tangent basis 선택과 독립적이다. Constraint도 같은 물리 공간으로 변환한 경우에만 spectrum의 좌표 불변성을 비교한다.

M6. 직접 적분 가능한 평판 oracle

평행 director와 affine 직사각형 면적 A에서는

[ C^{(0)}=\frac{\rho t A}{36} \begin{bmatrix}4&2&1&2\2&4&2&1\1&2&4&2\2&1&2&4\end{bmatrix}, \quad C^{(1)}=0,\quad C^{(2)}=\frac{t^2}{12}C^{(0)}. \tag{M10} ]

각 node는 mi=rho t A/4, hi=0, ji=rho t³ A/48이다. 이는 production quadrature를 호출하지 않고 다항 적분으로 얻는 oracle다. 곡면에서는 M5~M9와 independent volume integration을 사용한다. 두께 t 변경 시 stiffness와 rotary inertia는 다른 차수로 변하므로 전체 spectrum에 하나의 thickness power law를 강제하지 않는다.