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

9.2 KiB
Raw Blame History

질량과 공변 numerical stabilization

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. 비선형 objective potential과 고정 scale

기존 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을 생략하지 않는다. 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으로 출력하지 않는다.