Files
FESADev/docs/mitc4-shell-feast-eigen2/formulation-mass-stabilization.md

10 KiB
Raw Permalink Blame History

질량과 공변 numerical stabilization

현재 구현·검증은 기능 요약을 먼저 따른다. 아래 준비 상태·owner 인계와 옛 source 문서명은 작성 당시 이력이며, 원문은 Git 이력에서 확인한다.

Feature mitc4-shell-feast-eigen2, formulation r1. 기본 보고서의 적용·채택 경계를 따른다.

M1. 운동학과 physical stiffness

길이·질량·시간 차원은 L,M,T, 두께는 t, 시간은 tau다. 일정한 rho>0 [M/L³], centered constant t>0 [L], E>0 [M/(LT²)], -1<nu<1/2를 사용한다. 밀도 필수 조건은 eigen에만 적용한다. 기존 지원 geometry/director·section 조건을 유지한다. 원 node X_i, incident element가 공유하는 unit director d_i, 오른손 tangent frame (a_i,b_i,d_i)를 사용한다. Corner 순서는 (-1,-1),(1,-1),(1,1),(-1,1)이다.

[ 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}]. \tag{M1} ]

전역 q_i=[u_X,u_Y,u_Z,theta_X,theta_Y,theta_Z], u [L], theta [1]이며 양의 director 변화는 theta cross d다. S_i=-[d_i]_cross, P_i=I-d_i d_i^T=S_i^T S_i, G_i=[I,zS_i], H=[N_1G_1,...,N_4G_4]에서 v=H q_dot이다. 회전 physical 성분 alpha=a_i^T theta,beta=b_i^T theta, drilling gamma=d_i^T theta다. S_i theta=-b_i alpha+a_i beta, S_i d_i=0. 보간 director를 재정규화하지 않는다.

기존 shell 정식화의 geometry·MITC tying· 20↔24 변환을 그대로 사용한다. Engineering strain 순서 [epsilon11,epsilon22,gamma12,gamma13,gamma23]에 대해

[ 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}^0=\sum_{g\in2\times2\times2}\bar B_g^TC_5\bar B_gJ_gw_g, ] [ K_e=T_p^TK_{20}^0T_p+T_d^T(k_{d,e}I_4)T_d. \tag{M2} ]

T_p maps to 4×[uX,uY,uZ,alpha,beta]; T_d maps to four gamma values. Plane stress sigma33=0, 작은 변형·무응력 원형상이다. Parent momentum balance는 rho u_ddot=div sigma, 자유 traction은 0이다. MITC admissible variation의 약형은 integral rho delta u·u_ddot + integral delta e_bar^T C5 e_bar + delta gamma^T kd gamma=0. Lumping 후 반이산 잔차 r=M_L q_ddot+Kq, tangent는 각각 M_L,K이고 harmonic ansatz에서 K phi=lambda M_L phi, lambda=omega² [T^-2]다. Shear factor 5/6을 velocity/mass에 적용하지 않는다. 새 재료 history·감쇠·prestress· time integrator 또는 modal stress/reaction recovery는 없다.

M2. Consistent scalar moments와 적분

[ C^{(k)}{ij}=\int\rho N_iN_jz^kJ,d\xi d\eta d\zeta\quad(k=0,1,2),\qquad (M_C){ij}=\begin{bmatrix}C^{(0)}{ij}I&C^{(1)}{ij}S_j\ C^{(1)}{ij}S_i^T&C^{(2)}{ij}S_i^TS_j\end{bmatrix}. \tag{M3} ]

T_C=1/2 q_dot^T M_C q_dot. 각 C는 symmetric 4×4, 차원은 ML^k, mass block 차원은 M,ML,ML²다. M1의 J는 각 자연좌표 차수 ≤2이고 Ni Nj z^k J는 ≤(4,4,4)이므로 3×3×3 Gauss–Legendre로 정확 적분한다 (roundoff 제외). 축별 point [-sqrt(3/5),0,sqrt(3/5)], weight [5/9,8/9,5/9]다. 증가 xi→eta→zeta, node/component와 moment 합 순서를 고정한다. 모든 27점에서 기존 geometry predicate·finite·J>0를 검사한다. 이 검사를 전체 volume injectivity theorem으로 확대하지 않는다. Stiffness의 8점 적분은 유지한다. 곡면에서는 C1이 비영일 수 있고 C2=t² C0/12는 일반식이 아니다.

완성한 C와 M_C를 검증·보존하고 다음 row sum의 실제 원본으로 사용한다. Mass 원본을 lumped로 덮어쓰거나 독립 mass를 별도로 만들어 derived라 부르지 않는다.

M3. Derived nodal block lumping과 보존 범위

[ 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&h_iS_i\h_iS_i^T&j_iP_i\end{bmatrix},\quad M{L,e}=\operatorname{blockdiag}_i L_i. \tag{M4} ]

전체 24×24 row sum 또는 diagonal scaling이 아니라 scalar moment row sum 뒤 tensor 재구성이다. L_i=integral rho Ni Gi^T Gi dV이므로 node 간 coupling은 0, 같은 node의 6×6 coupling은 허용한다. Incident elements는 같은 d_i를 사용하고 stable source 순서로 block을 더한다.

[ q_i^TL_iq_i=m_i|u_i+(h_i/m_i)S_i\theta_i|^2+ (j_i-h_i^2/m_i)|S_i\theta_i|^2. \tag{M5} ]

양의 volume measure와 thickness에 의해 mi>0,ji>0,mi ji-hi²>0, rank=5, ker Li=span([0;di])다. 부호·finite·variance 검사를 clamp로 통과시키지 않는다. 거의 singular인 양의 mass를 rank cutoff로 제거하지 않고 conditioning 실패를 별도 진단한다.

보존하는 양은 sum mi=integral rho dV, sum hi=integral rho z dV, sum ji=integral rho z² dV, first spatial moment sum(mi Xi+hi di)=integral rho X dV다. Uniform translation energy는 consistent와 같다. ji Pi는 tangent rotary inertia이며 arbitrary rigid rotation의 orbital inertia tensor까지 consistent와 같게 만드는 계약은 아니다. 실제 차이는 vi=ui_dot+z Si theta_i_dot에 대해

[ \dot q^T(M_L-M_C)\dot q=\tfrac12\int\rho\sum_{ij}N_iN_j|v_i-v_j|^2dV\ge0. \tag{M6} ]

서로 다른 일반 속도장의 energy·spectrum 동일성은 요구하지 않는다. Affine rectangle 면적 A·평행 director의 독립 oracle은

[ 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=t^2C^0/12, \quad(m_i,h_i,j_i)=(\rho t A/4,0,\rho t^3 A/48). \tag{M7} ]

곡면 oracle는 독립 volume integration으로 C1와 변동 J까지 검사한다.

S1. 초기 trace 평균 선택과 공변성

R_i는 구속 전 초기 physical K20의 node i에 대한 2×2 tangent rotary principal block이다.

[ k_{ref,e}=\frac18\sum_{i=1}^4\operatorname{tr}R_i,\qquad k_{d,e}=10^{-3}k_{ref,e}>0. \tag{S1} ]

단위는 ML²/T²; 8개 rotary diagonal 전부를 평균한다. Translation diagonal·mass·density· trial tangent·BC·reference 오차는 입력이 아니다. Negative/nonfinite diagonal은 실패, exact zero는 평균에 포함, 8개 전부 zero면 positive inventory 실패다. 개별 Ri의 SPD를 새 geometry 조건으로 요구하지 않는다. Overflow 방지 계산은 dmax*((sum dr/dmax)/8)이며 고정 node1→4, alpha→beta 순서로 합한다. Underflow/nonpositive kd를 floor로 대체하지 않는다.

Proper global rotation Q와 재생성 tangent Ui'=Q Ui Oi에서 Ri'=Oi^T Ri Oi, tr Ri'=tr Ri. S_Q=blockdiag_i(Q,Q)이면 Kphys'=S_Q Kphys S_Q^T, Kdrill'=S_Q Kdrill S_Q^T, M'=S_Q M S_Q^T이다. Q=I인 독립 tangent basis 변경에도 성립한다. 새 geometry/frame에서 K20·kd·complete K를 모두 다시 계산해야 하며 이미 계산한 tensor만 회전해서는 이 선택을 검증할 수 없다. BC가 양쪽의 지원 component subset으로 같은 물리 free space를 표현할 때 constrained spectrum·static response의 공변성을 검증한다. 다른 BC를 같은 문제로 비교하거나 MPC를 추가하지 않는다.

물리 rigid field는 ui=c+Omega cross Xi, theta_i=Pi Omega이며 strain=gamma=0이다. 여섯 물리 rigid modes는 양의 mass를 갖는다. Pure numerical drill [0;di]는 mass=0, stiffness>0이다. Linear UR에 full spatial Omega를 넣는 field와 이 physical 대표값을 구분한다.

S2. 비선형 drilling 선택과 고정 scale

2026-10-20 현재 기본값은 초기 global 단위 director d_i0를 사용한 직접 스프링이다. gamma_i=d_i0·phi_i, Pi_d=0.5*kd*sum(gamma_i²), r_i=kd*gamma_i*d_i0, K_ii=kd*d_i0*d_i0^T이며 병진·절점 간 coupling은0이다. Drilling polar/swing/atan2를 평가하지 않고 유한 강체 회전의 인공 모멘트를 허용한다. S1의 초기 scale과 mass null은 유지한다. 원형상 선형/eigen block은 두 방식에서 같다. 다음 objective 유도는 선택 가능한 기존 상대 비틀림 방식에 한정한다.

기존 nonlinear 정식화 §10의 polar corotation, total nodal rotation, minimal swing 및 objective gamma_i=atan2(s_i,c_i)를 유지한다. S1을 초기 reference에서 한 번 계산하고 모든 trial·commit·increment·rollback에서 그 값을 고정한다. D kd=D² kd=0이다.

[ \Pi_d=\tfrac12 k_d\gamma^T\gamma,\quad r_d=k_dD_\gamma^T\gamma,\quad K_d=k_d\left(D_\gamma^TD_\gamma+\sum_i\gamma_iD^2\gamma_i\right). \tag{S2} ]

기존 방식에서는 polar/swing/atan2의 미분과 branch validation을 생략하지 않는다. 현행 구현의 polar1과 미분 정확도 실패를 보존하며 과거 polar12로 복원하지 않는다. 두 방식 모두 trial마다 scale을 바꾸지 않는다. 일반 nonlinear Hessian 전체의 SPD는 주장하지 않는다. 초기에도 Dgamma는 corotation coupling을 포함하므로 linear Td와 동일시하지 않는다. 기존 physical 선형극한, directional residual/Hessian, objective rigid orbit·state rollback을 검증한다. 새 eigen에는 이 nonlinear state를 전달하지 않는다. Physical strain/stress/reaction의 의미와 Newton·하중·commit 계약은 보존한다.

S3. 필수 partial-BC exact target

a=(0,1,0), b=(-4/5,0,3/5), d=(3/5,0,4/5), X=(-2a-b,2a-b,2a+b,-2a+b), t=1,E=12,nu=0,rho=1인 rectangle이다. 모든 DOF를 고정하되 node1 global URZ만 free다. Independent rotary block의 diagonal은 49/12,23/6, off-diagonal은 corner에 따라 ±1/8이다. S1에서

[ k_{ref}=95/24,\quad k_d=19/4800,\quad M_{ff}=3/50,\quad K_{ff}=10369/7500,\quad\lambda=10369/450. \tag{S3} ]

Global Z축 45° 회전 뒤에도 같은 BC를 유지하고 geometry/frame부터 재계산한다. URZ에 양의 mass가 있으므로 drill 좌표라고 삭제하면 안 된다. 이 값은 NR-MFE-001의 새 식 target이며 과거 pass 재사용이 아니다. NR가 독립 재유도한다. Static 대응 load는 node1의 tangent moment m=(5/3)M0 b, 다른 load=0이다. d·m=0,mZ=M0; theta1Z=M0/Kff, full residual r=Kq-F다. 회전 뒤 m'=Qm을 사용한다. Pure global eZ moment는 drilling 성분 때문에 기존 load 계약에 부합하지 않는다. Eigen의 constrained residual을 static reaction으로 출력하지 않는다.

2026-10-20 현재 이 static 사례의 직접 방식은 별도 인공 spatial 모멘트를 차감한 기존 gate를 통과하고, 기존 방식은 raw gate에서 계속 거부된다. S3의 K/M/lambda target이나 partial BC를 바꾼 결과가 아니며 Value 비교 예외를 정적 평형에 적용하지 않는다.