1. 서 론
2. 응용요소법의 정식화
3. 수치해석 및 결과 분석
3.1 강성 설정 전략별 벽체 해석
3.2 성분별 패치 테스트와 보정 계수 산정
3.3 보정 계수의 적용 범위 검증
3.4 격자 세분화에 따른 수렴 해석
3.5 스프링 강성 설정 지침
4. 결 론
1. 서 론
응용요소법(Applied Element Method, AEM)은 구조물을 다수의 강체 요소(rigid element)로 분할하고, 인접 요소를 그 경계면에 분포된 법선/전단 스프링으로 연결하여 연속체 거동을 근사하는 수치해석 기법이다(Meguro and Tagel-Din, 2000). 각 요소는 강체로 거동하므로 연속체의 변형과 응력은 전적으로 스프링을 통해 표현되며, 스프링 강성을 국소적으로 조정함으로써 균열의 발생·전파, 요소의 분리, 충돌·접촉, 그리고 붕괴까지를 비교적 단순한 모델로 추적할 수 있다(Tagel-Din and Meguro, 2000; Meguro and Tagel-Din, 2002). 이러한 장점으로 AEM은 폭발·충격·연쇄붕괴 해석을 위한 상용 소프트웨어 ELS(Extreme Loading for Structures; Applied Science International, 이하 ASI)에 구현되어 실무에서 널리 활용되고 있으며(Tagel-Din and Meguro, 2001; Grunwald et al., 2018; ASI, 2023), 철근콘크리트 구조의 비선형·반복하중 거동 해석에도 적용되어 왔다(Salem et al., 2011; Khalil, 2012). AEM은 본질적으로 Kawai(1978)의 강체–스프링 모델(rigid body spring model, RBSM) 및 lattice/spring 계열 방법(Bolander and Saito, 1998; Griffiths and Mustoe, 2001)과 같은 계보에 속한다.
AEM에서 연속체의 강성은 전적으로 스프링 강성으로 표현되므로, 스프링 강성 산정식의 선택은 FEM에서 재료 물성 입력만큼이나 해석 결과를 직접 지배한다. 그런데 이 선택에 대한 지침은 문헌과 매뉴얼 사이에서 일치하지 않는다. 초기 AEM 정식화(Meguro and Tagel-Din, 2000)는 전단 스프링 값으로 ks = GA / L를 제시하였으며, 이후 Tateo(2021)는 정사각형 RBSM 기본 셀의 변형에너지 등가로 ks = 2GA / L를 유도하였다. ELS 이론 매뉴얼(ASI, 2023)은 ks = 2GA / L를 기본값으로 채택하고, 얇은 휨 부재에 대해서는 휨–전단 보정(shear adjustment)을 추가로 권고한다. 즉 동일한 문제에 대해 사용자가 선택할 수 있는 강성 설정 전략이 여러 개 존재하며, 어느 것이 수렴된 기준해에 가까운 결과를 주는지에 대한 사용자 관점의 정량적 검증은 제한적이다.
AEM과 그 모태인 강체–스프링 계열 모델의 정확도 검증은 주로 두 방향에서 이루어져 왔다. 구조물 수준에서는 초기 정식화 연구가 캔틸레버 등 기본 예제로 이론해와의 일치를 보였고(Tagel-Din and Meguro, 2000), 폭발·연쇄붕괴 시나리오에서 FEM과 AEM의 예측을 모델 상세 수준별로 비교한 연구가 보고되었다(Grunwald et al., 2018). 정식화 수준에서는 Kawai(1978)의 RBSM 이후 스프링 강성의 유도가 요소 형상과 이산화에 의존한다는 점, 그리고 법선·전단 스프링만으로는 연속체의 Poisson 효과를 일관되게 표현하기 어렵다는 점이 lattice/스프링 계열 모델 일반의 제약으로 지적되어 왔다(Griffiths and Mustoe, 2001; Asahina et al., 2015). 그러나 이들 연구는 특정 정식화 하나의 검증이나 정식화 자체의 이론적 성질에 초점을 두었으며, 사용자가 실제로 선택해야 하는 복수의 강성 설정 — 초기 정식화, 매뉴얼 기본값과 보정, 보정 계수 — 을 동일 문제와 동일 기준해 위에서 비교하고 그 격자 의존성까지 추적한 연구는 찾기 어렵다.
이에 본 연구는 이 선택 문제를 사용자 관점에서 정량적으로 검토하고, 벽체류 문제에서 신뢰할 수 있는 강성 설정 방법을 도출한다. 횡력을 받는 2차원 벽체의 선형탄성 응답을 대상으로, 법선·전단 스프링 강성(kn, ks)을 산정하는 네 가지 전략을 동일 조건에서 비교한다. 네 가지 전략은 문헌·매뉴얼이 제시하는 ks의 세 가지 정식화(GA / L, 2GA / L, ELS shear-adjusted)와 성분별 기준문제인 패치 테스트(patch test)로 보정한 계수 α(법선)·β(전단)이며, 보정 성능의 상한을 보여주는 참조선으로 벽체 응답의 오차를 0으로 만드는 사후 계수 β를 함께 사용한다. 기준해로는 수렴 확인된 평면응력 FEM 해와 두 가지 정해(exact solution) 패치 테스트를 사용하고, 세장비 H / W = 2의 벽체에 대한 추가 해석으로 결과의 일반성을 확인한다.
논문의 구성으로 2장에서는 AEM 정식화를 설명하고, 3장에서는 강성 전략별 벽체 해석, 패치 테스트, 적용 범위 및 격자 수렴성을 검토한 뒤 실무적 설정 지침을 제시하며, 4장에서 결론을 정리한다.
2. 응용요소법의 정식화
AEM에서 구조물은 2차원의 경우 다수의 직사각형 강체 요소로 분할된다. 각 요소는 중심에 정의된 2개의 병진 자유도와 1개의 회전 자유도, 총 3 자유도를 가지며 요소 자체는 변형하지 않는 강체로 거동한다. 인접한 두 요소는 공유 경계면을 따라 분포된 법선·전단 스프링 쌍으로 연결되고, 법선 스프링은 접촉면에 수직, 전단 스프링은 평행하게 작용한다. Fig. 1(a)의 동일한 직사각형 요소 쌍에서 L은 요소 중심 간 거리, D는 접촉면 방향의 요소 치수, d는 개별 스프링 쌍의 분담폭이다. 면에 nₛ개의 스프링 쌍을 균등 배치하면 d = D / nₛ이고, 각 스프링 쌍의 대표면적은 A = d·t이다. 그림에 나타낸 수직 접촉면에서는 L = a, D = b이며, 정사각 요소에서는 L = D = h이다. Fig. 1(b)와 같이 두 요소를 잇는 스프링은 접촉점에서 양쪽 요소의 자유도를 연결하며, 요소 중심에서 접촉점까지의 레버암과 접촉면의 방향이 강성행렬을 결정한다.
한 스프링이 담당하는 법선 강성 kn과 전단 강성 ks는 대표 면적 A = d·t와 두 요소 중심을 잇는 스프링 길이 L로부터 식 (1)과 같이 산정된다. 여기서 E와 ν는 해석 대상 구조물 재료의 탄성계수와 Poisson비이고, G = E / [2(1 + ν)]는 그 전단탄성계수이다. 식 (1)의 전단 강성은 초기 AEM 정식화(Tagel-Din and Meguro, 2000)가 제시한 형태이며, 문헌 간 정의 차이와 그 영향은 3.1장에서 다룬다.
Fig. 1(b)와 같이 한 접촉점의 법선·전단 스프링 쌍은 양쪽 요소의 6개 자유도를 연결한다. 이들 자유도를 모은 벡터 q와, 접촉점의 위치 및 작용 방향으로부터 정해지는 법선·전단 변환 벡터 bn, bs를 식 (2)와 같이 정의한다.
식 (2)에서 bn과 bs는 각각 법선·전단 스프링의 기하 변환 벡터로, 동일한 형식에 작용 방향만 달리하여 구성된다. bn에는 접촉면의 단위 법선 n = (nx, ny)을, bs에는 그에 수직한 단위 접선 t = (–ny, nx)을 대입하며(t를 대입하면 식 (2)의 bs가 얻어진다), Fig. 1(b)의 θ는 요소 1의 중심과 접촉점을 잇는 선이 전역 x축과 이루는 각이다. (dx1, dy1)와 (dx2, dy2)는 Fig. 1(b)에 표시한 대로 각각 요소 1, 2의 중심에서 접촉점까지의 레버암이다. 각 스프링의 신장량 δn, δs는 해당 변환 벡터와 자유도 벡터의 곱으로 δn = bnᵀq, δs = bsᵀq와 같이 주어지며, 한 접촉면의 6 × 6 강성행렬과 전역 강성행렬, 지배방정식은 식 (3)과 같이 구성된다.
식 (3)에서 한 접촉면의 강성행렬 Ke는 면당 ns개 적분점 각각에 대해 법선 강성 kn에 대한 외적(outer product) 항 bnbnᵀ과 전단 강성 ks에 대한 외적 항 bsbsᵀ을 모두 합하여 구성된다(bbᵀ는 변환 벡터의 외적으로 얻어지는 계수(rank) 1의 행렬이다). 두 강성은 회전 자유도 항에서 레버암을 통해 함께 나타나며, 적분점마다 레버암이 달라 분포된 스프링들이 축력뿐 아니라 우력을 통해 휨모멘트를 전달한다. 이것이 단일 스프링이 아니라 면당 다수의 분포 스프링을 사용하는 이유이다. 모든 접촉면의 Ke를 합하여 전역 강성행렬 K를 얻고, 경계조건을 부여한 뒤 식 (3)의 지배방정식을 풀어 요소 변위를 구한다. 이와 같이 외적 합으로 표현한 Ke는 Meguro and Tagel-Din(2000)이 제시한 삼각함수로 표현된 닫힌형 요소강성행렬과 수학적으로 동등하며, 본 연구는 면당 분포 스프링으로의 일반화가 용이한 외적 합 형태를 채택하였다. FEM에서 절점이 요소로 연결되는 것과 달리 AEM에서는 요소가 스프링으로 연결되며, 미지수가 절점이 아니라 요소에 정의되므로 전역 자유도 수는 요소 수에만 의존한다.
AEM에서는 자유도가 요소 중심에 정의되고 강성행렬도 요소 중심에 가한 단위 변위로부터 구성되므로, 본 연구에서는 경계 요소의 중심 자유도(u1, u2, u3)를 직접 구속하여 지지조건을 부여하였다(Tagel-Din and Meguro, 2000). 면당 스프링 수 ns는 분포 강성과 휨모멘트 전달의 수치적 표현 정밀도를 결정하는 이산화 파라미터이며, 스프링 수를 늘려도 전역 자유도는 변하지 않아 계산비용 증가가 작다. 비선형 해석에서는 각 스프링 강성을 재료 모델에 따라 갱신하고 인장강도 초과 시 스프링을 제거하여 균열·분리·붕괴를 모사하지만, 본 연구는 선형탄성 거동에 한정한다.
3. 수치해석 및 결과 분석
3장에서는 2장의 AEM 정식화를 바탕으로 스프링 강성 설정 전략을 정의하고, 이를 검증하기 위한 네 종류의 수치해석을 수행한다. 모든 해석은 횡력을 받는 캔틸레버형 2차원 철근콘크리트 벽체 패널, 즉 폭 W = 3.0 m, 높이 H = 3.0 m, 두께 t = 0.20 m의 정사각 패널을 기본 해석 대상으로 하며, 콘크리트 Ec = 30 GPa, ν = 0.20, 철근 탄성계수 Es = 200 GPa, 분산 보강비 ρx = ρy = 0.002로 모델링하고 절마다 경계조건과 하중만 달리한다. AEM 해석은 모든 경우 면당 5개의 스프링 쌍을 사용하고, FEM·AEM 공통으로 강성행렬을 희소행렬 형식으로 조립하여 Jacobi 전처리 공액구배법(PCG; Saad, 2003)으로 푼다. 반복 수렴의 허용오차(tolerance)는 상대 잔차 10–10로 설정하였다. 오차는 (해석값 − 기준값) / 기준값 × 100(%)로 정의한다. 요소 크기 h는 정사각 요소 한 변의 길이를 뜻하며, 오차가 h에 비례하여 감소하는 1차 수렴 거동을 O(h)로 표기한다.
3.1 강성 설정 전략별 벽체 해석
본 절에서는 횡력을 받는 정사각 벽체 한 문제를 대상으로, 문헌·매뉴얼의 세 가지 강성 설정 전략(S1~S3)이 정점 변위의 정확도에 미치는 영향을 다섯 가지 격자에 걸쳐 비교한다. 이를 위해 먼저 본 연구가 다루는 전략 전체를 정의하고, 벽체의 FEM 기준해를 수렴 확인한 뒤, 전략별 AEM 해석 결과를 분석한다. 패치 보정 전략(S4)은 보정 계수를 3.2장에서 산정한 뒤 3.3장에서 평가한다.
2장의 일반식에서 법선 강성은 문헌과 매뉴얼이 동일하게 kn = E·A / L를 제시하는 반면, 전단 강성은 GA / L과 2GA / L로 2배 갈린다. 분산철근은 콘크리트 총단면에 방향별 철근 강성을 가산하여 반영하였다. 즉, E = Ec + ρᵢEₛ = 30 + 0.002(200) = 30.40 GPa이고 G = Ec / [2(1 + ν)] = 12.50 GPa이다. 본 연구는 법선 스프링을 kn = α·E·A / L, 전단 스프링을 ks = β·2G·A / L 형태로 두고(식 (4)), 보정 계수 α(법선)와 β(전단)로 네 가지 설정 전략을 Table 1과 같이 통합 표현한다. (S1) β = 0.5는 초기 정식화의 GA / L에 해당한다(Meguro and Tagel-Din, 2000). (S2) β = 1.0은 ELS 기본값 2GA / L이다(ASI, 2023). (S3) ELS shear-adjusted는 얇은 휨 부재를 겨냥해 2GA / L에 Table 1의 휨–전단 보정계수(flexural correction)를 곱하는 매뉴얼 권고 전략으로(ASI, 2023), 보정계수는 요소 스트립의 전단 강성과 휨 강성을 직렬 결합해 유도된 1 / [1 + 2G·L2 / (E·D2)]이다. 여기서 L은 식 (1)의 스프링 길이, D는 요소 스트립의 휨 깊이(경계면 방향의 요소 치수)이다. 정사각 요소에서는 L = D가 되어 보정계수가 격자와 무관한 상수 (1 + ν) / (2 + ν)(ν = 0.2에서 약 0.55)로 고정되므로 유효 전단 강성은 약 1.09·GA / L이 되며, L ≪ D인 스트립에서만 1에 수렴해 2GA / L로 돌아간다. (S4) 패치 보정은 각 스프링 성분을 고립하는 기준문제에서 α, β를 3.2장에서 격자별로 산정해 적용하는 전략이다. 아울러 벽체 정점 응답 자체를 기준해에 맞추는 사후(post hoc) 계수 β를 함께 산정하는데, 이는 실무에서 쓸 수 있는 전략이 아니라 보정이 도달할 수 있는 상한 성능을 보여주는 참조선이다.
Table 1
Spring-stiffness setting strategies compared in this study
| Strategy | Shear spring ks | Normal spring kn | Idea / Source |
| (S1) AEM-2000 basic | GA / L | EA / L | Meguro and Tagel-Din (2000) |
| (S2) ELS basic | 2GA / L | EA / L | ASI (2023) default |
| (S3) ELS shear-adjusted | 2GA / L × flexural correction | EA / L | ASI (2023) guideline for thin members |
| (S4) Patch-calibrated | β·2GA / L | α·EA / L | This study (exact patches) |
전략 비교의 대상 문제는 하단을 완전 고정하고 우측 변에 합력 60 kN의 등분포 수평 표면력(qx = 100 kPa)을 가한 횡하중 벽체이며, 주 응답량은 Fig. 2에 빨간 점으로 표시한 우측 상단 모서리의 수평변위 ux(이하 정점 변위)이다. 기준해 산정을 위한 FEM 해석은 벽체를 4절점 평면응력 사각형 요소로 이산화하고, 분산 보강은 앞서와 같은 혼합법칙으로 재료 강성행렬의 두 축방향 대각 성분에 반영하였다. 격자는 10 × 10에서 50 × 50까지 5종(자유도 242~5,202)을 사용하고, 우측 변의 등분포 표면력은 절점 등가하중으로 치환하였다. Table 2는 격자별 정점 변위와 50 × 50 해 대비 오차를 정리한 것이다. FEM 오차는 단조 감소하여 40 × 40에서 50 × 50 격자 대비 약 0.033 %에 불과하므로, 50 × 50 해(정점 변위 0.03331 mm)를 수렴된 단일 기준해(이하 FEM50)로 채택하였다. 모든 격자에서 작용 하중과 반력이 기계정밀도로 평형을 이루어 모델링의 일관성을 확인하였다.
Table 2
FEM reference: wall mesh convergence (vs FEM50)
| FEM mesh | DOF | Wall tip ux (mm) | Relative error of ux (%) |
| 10 × 10 | 242 | 0.032967 | –1.019 |
| 20 × 20 | 882 | 0.033218 | –0.265 |
| 30 × 30 | 1,922 | 0.033273 | –0.098 |
| 40 × 40 | 3,362 | 0.033295 | –0.033 |
| 50 × 50 | 5,202 | 0.033306 | 0.000 |
AEM 해석은 벽체를 10 × 10, 15 × 15, 20 × 20, 30 × 30, 40 × 40의 정사각 강체 요소로 분할하고(자유도 300~4,800), 하단 요소열의 중심 자유도를 완전 구속하고 우측 요소열의 요소 중심에 표면력 등가 수평하중을 가하여 수행하였다. 이 조건에서 전략 S1~S3 각각에 대해 격자별 정점 변위를 계산하여 FEM50과 비교하였다. Table 3은 세 전략의 격자별 정점 변위 오차를 정리한 것이고, Fig. 3(a)는 이를 격자 수에 대해 나타낸 것이다. 30 × 30 격자에서 전단 정식화 선택만으로 정점 변위 오차가 –8.6 %(2GA / L)에서 +30.1 %(GA / L)까지 변한다. GA / L은 전단 강성을 절반으로 두어 변위를 크게 과대평가하고, ELS shear-adjusted(+23.7 %)도 같은 방향으로 크게 빗나간다. 즉 스프링 강성 설정은 이 문제에서 결과를 지배하는 민감 입력이다.
더 중요한 것은 Fig. 3(a)와 Table 3에 나타난 격자 의존성이다. GA / L 계열 두 전략(S1·S3)의 오차는 10 × 10에서 +6.5~+12.1 %로 오히려 작아 보이지만, 격자를 세분화할수록 +26~+32 %까지 단조 증가한다. 즉 이들은 FEM 해가 아닌 다른(과유연한) 극한으로 수렴하며, 거친 격자에서의 근사한 일치는 이산화 오차(과강성)와 정식화 오차(과유연)가 상쇄된 우연의 결과다. 특히 매뉴얼이 얇은 부재에 권고하는 shear-adjusted(S3)는 10 × 10에서 +6.5 %로 세 전략 중 가장 정확해 보이나 40 × 40에서는 +25.9 %로 반전된다. 격자 수렴 확인 없이 거친 격자의 일치만으로 전략을 선택하면 오판하게 됨을 보여주는 결과다. 반면 2GA / L(S2)만이 세분화와 함께 FEM 기준해로 접근한다.
Table 3
Wall tip error (%) of each strategy versus mesh (vs FEM50)
이 격자 의존성이 특정 형상에 국한되지 않음을 확인하기 위해, 폭 W = 3.0 m, 높이 H = 6.0 m의 세장 벽체(H / W = 2)에 대해 동일한 비교를 추가로 수행하였다. 정사각 요소가 유지되도록 격자는 10 × 20, 20 × 40, 30 × 60을 사용하고, 하중은 우측 변의 동일 표면력(합력 120 kN), 기준해는 같은 방식으로 수렴을 확인한 FEM 50 × 100 해(정점 변위 0.2994 mm, 40 × 80 격자 대비 차이 0.02 %)이다. Table 4와 Fig. 3(b)는 세장 벽체에 대한 같은 비교의 결과이다. GA / L 계열은 가장 거친 격자에서 +1.5~+4.3 %로 오히려 작아 보이지만 세분화와 함께 +10~+13 %까지 증가하고, 2GA / L은 –12.6 %에서 –4.7 %로 단조 수렴한다. 즉 위의 경향은 세장비가 달라져도 유지된다.
Table 4
Wall tip error (%) of each strategy for the slender wall (H / W = 2, vs FEM 50 × 100)
마지막으로 법선·전단 성분의 상대 민감도를 보기 위해 법선 강성을 kn = EA / L(α = 1)로 고정하고 β만 변화시킨 결과가 Fig. 4(a), 반대로 전단 강성을 ks = 2GA / L(β = 1)로 고정하고 α만 변화시킨 결과가 Fig. 4(b)이다. 두 경우 모두 세 격자에서 정점 변위 오차가 계수에 대해 단조 감소하며, Fig. 4(a)에서 오차는 β = 0.5(GA / L)의 큰 양의 오차에서 β = 1.5의 음의 오차까지 단조 감소하며, 기준선으로 표시한 β = 1.0(2GA / L)에서도 세 격자 모두 음의 오차를 보인다. 오차를 0으로 만드는 β는 격자마다 달라 하나의 β로 모든 격자를 맞출 수 없음을 보여준다. 두 패널을 비교하면 같은 폭으로 잘못 설정했을 때 응답을 흔드는 정도는 법선이 전단보다 오히려 약간 크지만, 법선은 문헌 간 제시식이 일치해 사용자가 선택할 여지가 없는 반면 전단은 2배 갈리므로, 실질적인 선택의 문제는 전단에 있다. 요컨대 Fig. 4는 두 가지를 보여준다. 첫째, α, β의 선택에 따라 오차가 수십 %까지 달라지므로 스프링 강성 계수는 해석 결과를 지배하는 민감 입력이다. 둘째, 오차를 0으로 만드는 계수는 격자마다 이동하므로, 특정 격자에 맞춘 보정 없이 일반적으로 쓸 수 있는 선택은 α = β = 1이며 이때의 잔여 오차는 3.4장과 같이 격자 세분화로 관리한다.
3.2 성분별 패치 테스트와 보정 계수 산정
패치 테스트(patch test)는 이산화가 상수 변형률(응력) 상태를 정확히 재현하는지를 통해 정식화의 일관성과 수렴의 필요조건을 검증하는 고전적 시험이다(Bazeley et al., 1965; Taylor et al., 1986). 본 절에서는 같은 취지로 3.1장에서 확인한 오차의 근원을 성분별로 파악하기 위해, 동일 패널에서 법선 스프링만 활성화되는 축방향 패치와 전단 스프링만 활성화되는 순수전단 패치의 두 정해 기준문제를 다룬다. 각 패치로 보정 없는 상태의 오차를 정량화하고, 격자별로 오차를 0으로 만드는 보정 계수(이하 영오차 계수, zero-error factor) α*, β*를 산정한다.
두 패치 문제의 구성은 Fig. 5와 같다. Fig. 5(a)의 축방향 패치는 좌측 변 요소열의 수평 자유도를 롤러로 구속하고 좌하단 요소의 연직 자유도를 추가 구속한 뒤, 우측 변에 qx = 100 kPa의 등분포 표면력을 가한 것이다. 패널 전체가 상수 단축응력 상태가 되어 변형이 법선 스프링에만 집중된다. 이때 우측 변 수평변위의 정해는 u = qx·W / Ē이며, 여기서 Ē는 분산철근과 평면응력 Poisson 결합을 포함한 축방향 유효계수로, Ē = D11 – D122 / D22 = 30.42 GPa이며, u = qₓW / Ē = 9.8633 μm이다. Fig. 5(b)의 순수전단 패치는 네 변에 자기평형을 이루는 전단 표면력 τ = 100 kPa을 가하고 좌하단 모서리 요소의 세 자유도만 구속한 것으로, 상수 전단응력 상태가 되어 전단 스프링만 활성화되며, 변형에너지의 정해는 U = τ2·W·H·t / (2G) = 0.7200 J이다. 축방향 패치는 우측 변 수평변위가 위치에 무관한 단일 값이어서 변위를 기준으로 삼은 반면, 순수전단 패치는 상수 전단응력에 대응하는 변위장이 강체 회전의 구속 방식에 따라 달라지는 관례 의존량이므로 회전에 불변인 변형에너지를 기준으로 삼았다. 두 상태 모두 상수변형률이므로 4절점 FEM은 가장 거친 격자에서도 정해를 기계정밀도로 재현함을 확인하였으며, 이 정해를 성분별 기준으로 사용한다.
두 패치를 10 × 10~40 × 40의 정사각 AEM 격자로 해석하여, 보정이 없는 상태(α = β = 1)의 오차와 각 격자에서 오차를 0으로 만드는 성분별 영오차 계수 α*, β*를 구하였다. α*는 축방향 패치의 변위 오차가, β*는 순수전단 패치의 에너지 오차가 부호를 바꾸는 구간을 이분법으로 좁혀 ±0.001 이내로 결정하였다. Table 5가 그 결과이다.
Table 5
Normal- and shear-patch errors and zero-error factors on square meshes
Table 5는 두 가지를 보여준다. 첫째, 보정이 없으면 두 스프링 모두 유한 격자에서 과강성이다. 법선의 오차는 10 × 10의 –9.95 %에서 40 × 40의 –2.45 %로, 전단(2GA / L)의 오차는 –5.34 %에서 –1.29 %로 줄어들어, 둘 다 격자 세분화와 함께 소멸하는 이산화 오차임을 알 수 있다. 전단 오차가 이렇게 0으로 수렴한다는 것은 2GA / L이 점근적으로 옳은 전단 강성이라는 뜻이고, 뒤집어 말하면 그 절반인 GA / L은 세분화로도 사라지지 않는 정식화 오류를 안고 있다는 뜻이다. 3.1장에서 GA / L 계열이 잘못된 극한으로 수렴한 근원이 여기서 확인된다. 둘째, 영오차 계수는 법선 α*가 0.900에서 0.976으로, 전단 β*가 0.946에서 0.987로 격자 세분화와 함께 1.0에 접근하고, 이 값들에서 패치 오차는 1 % 미만이 된다. 즉 보정 계수는 본질적으로 이산화 오차의 보상이어서 격자가 가늘어질수록 보정할 것이 없어지며, 특정 격자에서 얻은 값은 그 격자에서만 유효하다.
이 이산화 오차가 소멸하는 속도는 Fig. 6으로 확인된다. Fig. 6은 α = β = 1의 두 패치 오차를 자유도 수에 대해 log–log로 나타낸 것으로, 정사각 격자에서 자유도가 3n2(n은 변당 요소 수)으로 h–2에 비례하므로 기울기 –1 / 2의 직선이 O(h)에 해당한다. 두 패치의 오차가 모두 이 기울기를 따르므로, 법선·전단의 과강성은 요소 크기에 비례해 1차 속도로 소멸하는 O(h) 오차임이 확인된다.
3.3 보정 계수의 적용 범위 검증
본 절에서는 3.2장에서 산정한 보정 계수의 적용 한계를 검증하기 위해 두 가지 수치해석 예제를 다룬다. 요소 종횡비를 변화시킨 순수전단·축방향 패치와, 패치 보정 계수를 그대로 적용한 벽체 횡하중 해석이다.
요소 종횡비의 영향은 다음 해석으로 검토하였다. 동일 패널을 24 × 16, 40 × 20, 30 × 10으로 분할하면 요소 종횡비 dx / dy가 각각 0.67, 0.50, 0.33이 된다. 이들 격자와 정사각 기준인 20 × 20(dx / dy = 1.00)에 대해 순수전단 패치와 축방향 패치를 모두 해석하여, 보정이 없는 상태(α = β = 1)의 오차와 성분별 영오차 계수를 구하였다. 종횡비가 1에서 먼 경우 β*가 통상 범위를 벗어날 수 있어 β*의 탐색 범위는 0.5–16으로 확장하였다. Table 6이 그 결과이고, Fig. 7은 이를 종횡비에 대해 나타낸 것이다. 정사각 요소에서 –2.6 %였던 순수전단 패치의 β = 1 오차는 종횡비가 1에서 멀어질수록 +19.8 %, +57.4 %, +114.3 %로 급증하고, 이를 상쇄하는 β*도 0.97에서 1.24, 2.05, 9.42로 커져 물리적 보정의 범위를 완전히 벗어난다. 반면 같은 요소에서 축방향 패치의 α = 1 오차는 –5.0 %에서 –2.4~–4.1 %로 종횡비에 거의 반응하지 않고, α*도 0.950–0.976으로 1.0 부근을 유지한다. 즉 3.2장의 깨끗한 O(h) 거동은 정사각 요소에 한정되며, AEM의 전단 표현은 요소 종횡비에 매우 취약한 반면 법선 표현은 강건하다. 이러한 법선 표현의 강건성은 문헌과 매뉴얼이 법선 강성을 이견 없이 EA / L로 제시해 온 것과 부합하며, 강성 설정 논쟁이 전단에 집중되어 온 이유를 보여준다.
Table 6
Shear- and normal-patch errors and zero-error factors versus element aspect ratio
Table 7은 3.1장과 동일한 벽체 횡하중 문제에 보정 없음(ELS basic), 패치 보정, 정점 영오차 β의 세 가지 강성 설정을 적용했을 때의 계수와 정점 변위 오차를 정리한 것이다. 패치 보정에는 Table 5에서 격자별로 산정한 영오차 계수 (α*, β*)를 그대로 사용하여 성분별 보정이 다른 문제로 전이되는지를 검토하였고, 정점 영오차(zero tip error) β는 법선 강성을 kn = EA / L(α = 1)로 고정한 채 벽체 정점 변위 오차를 0으로 만드는 사후 계수로서 동일한 이분법으로 구하였다. ELS basic의 오차는 Table 3의 S2와 같은 값이고, 영오차 설정의 오차는 정의상 0이라 생략하였다. 패치 보정 전략(S4)도 정점 변위 오차가 10 × 10 – 14.6 %, 40 × 40 – 5.1 %로 여전히 크다(보정 없는 ELS basic은 각각 –21.5 %, –7.0 %). 정점 영오차 β는 0.61–0.85로 순수전단 패치의 β*(0.95–0.99)와 뚜렷이 다르다. 짧은 캔틸레버의 지배 모드인 휨은 축·전단 패치 어디서도 직접 보정되지 않기 때문이다. 보정 계수는 산정에 사용한 상수응력 모드에서는 정의상 정확하지만 다른 변형 모드로는 전이되지 않으므로, 영오차 보정 계수는 재료 상수가 아니라 하중·변형 모드에 의존하는 양이다.
Table 7
Wall lateral-load response under three factor settings (vs FEM50): uncalibrated ELS basic, patch-calibrated, and zero-error
3.4 격자 세분화에 따른 수렴 해석
본 절에서는 보정 없는 ELS basic(β = 1) 설정 하나만을 대상으로, 벽체 횡하중 문제의 격자를 10 × 10에서 160 × 160까지 세분화하는 수렴 해석을 다룬다. 이로부터 격자 세분화 경로로 도달할 수 있는 정확도와 남는 한계를 정량화한다.
해석은 벽체 격자를 10 × 10, 15 × 15, 20 × 20, 30 × 30, 40 × 40, 60 × 60, 80 × 80, 120 × 120, 160 × 160(요소 25,600개, 자유도 76,800)의 9단계로 세분화하며, 각 격자에서 α = β = 1 설정으로 정점 변위를 계산하는 방식으로 수행하였다. 정확도의 기준해는 어디까지나 FEM50이다. 다만 오차를 격자 세분화로 사라지는 이산화 성분과 사라지지 않는 모델 성분으로 분해하기 위해, AEM 해 자체가 수렴하는 극한 u∞를 진단량으로 함께 추정하였다. u∞는 격자비 2의 세 격자(40 × 40, 80 × 80, 160 × 160)의 해에, 수렴 극한과 수렴차수를 근사적으로 구하는 표준 기법인 Richardson 외삽을 적용하여 다음과 같이 얻었다(Richardson, 1911; Roache, 1994).
Fig. 8(a)의 0.69와 0.99는 각각 FEM50과 u∞를 기준으로 한 오차를 e = Chp로 두고, 모든 격자의 (log h, log e) 자료를 최소자승 적합하여 얻은 수렴차수 p이다. Fig. 8(b)는 3.3장의 정점 영오차 β도 격자별로 함께 산정하여 격자 세분화에 따른 이동을 나타내었다. 정점 변위 오차는 FEM50 대비 –21.5 %(10 × 10)에서 –3.2 %(160 × 160)로 단조 감소하고, 정점 영오차 β도 0.61에서 0.93(160 × 160)으로 1.0을 향해 표류하되 도달하지는 않는다. 그러나 개선 폭은 갈수록 줄어들어, 요소 크기를 절반으로 줄이는 80 × 80에서 160 × 160의 세분화가 오차를 1.3 %p 줄이는 데 그친다. FEM 기준의 겉보기 수렴률도 약 h0.69로 1차에 못 미치지만, AEM 해 자체의 Richardson 외삽 극한(u∞)을 기준으로 하면 수렴률은 약 h0.99로 패치 테스트에서 확인한 O(h) 오차와 일치한다. 대신 이 극한은 FEM 연속체 해보다 약 1.9 % 작은(과강성) 값으로, 법선·전단 스프링만으로는 Poisson 결합을 표현할 수 없는 격자–스프링 계열 모델의 알려진 한계(Asahina et al., 2015)에 대응하며, 이 잔여 오차는 격자 세분화로도 사라지지 않는다. 즉 충분한 세분화는 보정 없이도 해를 연속체 수준에 접근시키지만, 보정 계수는 거친 격자에서 일치하는 모드의 정확도를 높이는 가속 인자일 뿐 격자 세분화를 대체하지 못한다.
덧붙여 이 수렴 거동은 일반적 보장이 아니라 본 예제에 대한 경험적 관찰이다. FEM은 변분 정식화 위에서 수렴이 수학적으로 증명되고 상수변형률 패치 테스트를 통과하지만(본 예제의 FEM 수렴률 약 h1.7; Zienkiewicz et al., 2013), AEM에는 그러한 일반 수렴 이론이 없다. 따라서 새로운 문제 유형에 AEM을 적용할 때는 격자 수렴 확인이 항상 뒤따라야 한다.
3.5 스프링 강성 설정 지침
이상의 결과를 종합하면, 횡력을 받는 2차원 벽체의 AEM 해석에서 스프링 강성 설정과 격자 운용에 대해 다음의 지침을 제시할 수 있다.
(1) 전단 스프링 강성은 ks = 2GA / L(β = 1)로 설정한다. 초기 문헌의 GA / L은 강성을 절반으로 두는 잘못된 정식화이며, 이를 쓰는 전략은 격자를 세분화할수록 FEM 해에서 멀어진다.
(2) ELS shear-adjusted와 같은 휨–전단 보정은 벽체류 문제에는 적용하지 않는다. 거친 격자에서 정확해 보이는 것은 이산화 오차와 정식화 오차가 상쇄된 착시이며, 격자 수렴 확인 없이 거친 격자 결과로 전략을 선택해서는 안 된다.
(3) 요소는 가능한 한 정사각형으로 유지한다. AEM 전단 표현은 요소의 종횡비에 민감하며, 종횡비 0.5인 요소에서는 β*가 약 2.1, 요소 종횡비 0.33에서는 약 9.4에 달해 물리적 보정 범위를 벗어난다.
(4) 정확도는 보정 계수가 아니라 격자 세분화로 확보하고, 반드시 두 개 이상의 격자로 수렴을 확인한다. 본 벽체 예제에서 이산화 오차는 요소 크기에 비례하는 대략 –200 % × h / W여서 목표 오차에 필요한 분할 수를 가늠하는 출발점이 될 수 있다. 다만 이 비례상수는 구조물의 형태와 변형 모드에 따라 달라질 수 있으므로 다른 문제에 일반화할 수 없으며, 격자 수렴 확인을 대체하지도 않는다.
(5) 패치 보정 계수는 하중이 보정에 사용한 상수응력 모드와 일치할 때만 거친 격자의 가속 수단으로 제한적으로 사용한다. 휨 지배 문제로는 전이되지 않아 5–15 %의 오차가 남는다.
(6) 위 설정을 모두 따라도 Poisson 결합 부재로 인한 약 2 %의 과강성 오차는 남는다. 이보다 높은 정확도가 요구되는 선형탄성 문제라면 AEM보다 FEM이 적합하다.
4. 결 론
본 연구는 횡력을 받는 2차원 벽체의 AEM 해석 정확도가 스프링 강성 설정에 의해 어떻게 지배되는지를 문헌·매뉴얼의 세 가지 정식화와 패치 보정 전략의 통제 비교로 정량화하고, 수렴된 FEM 기준해에 근접하기 위한 설정 방법을 도출하였다. 이를 기반으로 다양한 수치해석을 수행하였으며, 다음과 같은 결론을 얻었다.
1) 동일 격자에서 전단 정식화 선택만으로 정점 변위 오차가 –8.6 %(2GA / L)에서 +30.1 %(GA / L)까지 변하였다. GA / L 계열(초기 정식화, ELS shear-adjusted)은 격자 세분화와 함께 오차가 +26~+32 %까지 증가하며 잘못된 극한으로 수렴하였고, 거친 격자에서의 일치는 이산화 오차와 정식화 오차가 상쇄된 결과였다. 정해 패치 테스트는 2GA / L이 점근적으로 옳은 전단 강성임을 확인해 주며, 이 경향은 세장비 H / W = 2의 벽체에서도 동일하게 확인되었다.
2) 법선 kn = E·A / L와 전단 ks = 2G·A / L 모두 유한 격자에서 O(h) 과강성 오차를 가지며, 패치 보정 계수 α·β는 격자 세분화와 함께 1.0으로 수렴하지만 휨 지배 벽체로 전이되지 않았고(잔여 5–15 %), 보정에 사용한 변형 모드에서만 유효하였다. 즉 영오차 보정 계수는 재료 상수가 아니라 하중·변형 모드, 격자, 요소 형상에 의존한다.
3) AEM 모델링에서 기본적으로 신뢰할 수 있는 방법은 i) 정사각 요소 사용; ii) 2GA / L(β = 1) 전단 스프링 강성 사용; iii) 격자 세분화 수행이다. 이 경로에서 AEM 해는 자체 극한으로 약 O(h)로 수렴하나, 그 극한은 Poisson 결합 부재로 인해 연속체 해보다 약 2 % 강성측에 있으며(FEM 기준 겉보기 수렴률 약 h0.69), FEM과 같은 일반 수렴 이론이나 패치 테스트 통과는 보장되지 않는다.
4) 따라서 벽체류 문제의 해석에서는 전단 강성을 2GA / L로 설정하고, 휨–전단 보정은 벽체류 문제에 사용하지 않으며, 요소는 정사각을 유지하고, 정확도는 격자 수렴 확인으로 관리하는 것이 타당하다. 보정 계수는 하중 모드가 명확한 경우의 거친 격자 가속 수단으로만 사용한다.
첨언으로 본 연구는 선형탄성·정적 범위에 한정되며, 비선형·동적 확장과 Poisson 결합의 일관된 정식화는 후속 연구로 남긴다.
