OpenSees 쉘·솔리드·지반 요소: 같은 요소 10개가 77% 뻣뻣한 답부터 99% 무른 답까지 내놓는 이유

같은 캔틸레버, 같은 요소 10개, 같은 하중, 그런데 답은 77% 뻣뻣한 것부터 99% 무른 것까지 다섯 가지입니다. 바뀐 것은 element 뒤의 단어 하나뿐입니다. 이 시리즈의 지금까지 모든 편에서 요소는 보였고, 보 요소는 너그럽습니다. 2절점 elasticBeamColumn은 등단면 부재를 요소 하나로 정확히 풉니다. 연속체 요소는 너그럽지 않습니다. 메시가 모델의 일부이고, 요소 정식화가 모델의 일부이며, 둘 다 틀렸을 때 스스로 알려주지 않습니다.

이 글은 1차원 요소를 떠날 때 무엇이 달라지는지를 다룹니다. 쉘, 평면요소, 솔리드, 지반. 모든 수치는 OpenSeesPy 3.8.0에서 산출했고 닫힌해 또는 공표된 벤치마크와 대조했습니다. 이 글은 11부입니다. 1–10부에서 객체 모델, 탄성 골조, 파이버 단면, 동적해석, 소스 읽기, 철근콘크리트, 지진기록, 3차원 건물, 수렴 실패, Tcl 읽기를 다뤘습니다.

1. 달라지는 세 가지

수치에 들어가기 전에, 골조 모델에서 가져온 습관을 깨는 차이 셋입니다.

골조 요소연속체 요소

재료uniaxialMaterial, 응력 성분 하나nDMaterial, 응력 텐서 전체

절점당 자유도2D 3개, 3D 6개평면 2, 솔리드 3, 쉘 6

요소 하나대개 정확절대 정확하지 않음. 메시가 입력이다

헷갈리는 오류를 만드는 건 ndf 줄입니다. 평면응력 모델은 model('basic','-ndm',2,'-ndf',2), 솔리드 모델은 -ndm 3 -ndf 3, 쉘 모델은 쉘 절점이 회전을 가지므로 -ndm 3 -ndf 6입니다. 한 모델에서 쉘과 솔리드를 섞으려면 절점 정의 사이에서 ndf를 바꿔야 합니다. 합법이고, 이 영역에서 나오는 "invalid number of arguments" 메시지의 대부분이 여기서 옵니다.

2. 쉘: 공식 벤치마크에 요소 1000개가 필요하다

핀치드 실린더는 표준 쉘 문제 중 하나이고, OpenSees는 이것을 EXAMPLES/verification/PinchedCylinder.tcl로 배포합니다. 반경 300, 두께 3인 실린더를 마주 보는 두 점하중으로 누르고, 공표된 급수해가 하중점 처짐을 w = −164.24 P / (E t)로 줍니다. 이 스크립트를 파이썬으로 옮기고 메시 스윕을 확장한 결과입니다.

메시요소 수ShellMITC4ShellDKGQ

4 × 41665.01%36.09%

8 × 86427.17%4.97%

16 × 162567.49%1.55%

32 × 321,0241.23%1.05%

48 × 482,3040.11%0.69%

64 × 644,0960.64%0.51%

96 × 969,2161.09%0.35%

이 표에서 가져갈 것이 셋입니다.

4 × 4 쉘 메시는 성긴 메시가 아니라 틀린 메시입니다. 65%는 무엇의 근사도 아닙니다. 4분원에 요소 16개는 화면에서 그럴듯해 보이고, 그게 함정입니다.

ShellDKGQ는 요소 64개로 5%에 도달하는데 ShellMITC4는 400개쯤 필요합니다. 휨을 받는 얇은 쉘에서 이산 Kirchhoff 요소는 메시 기준 한 자릿수의 값어치가 있습니다. ShellNLDKGQ는 여기서 ShellDKGQ와 모든 자리가 같았습니다. 선형 해석에서는 당연합니다. 기하비선형 정식화가 붙었을 뿐 같은 요소입니다.

ShellMITC4가 가장 좋아 보이는 메시는 수렴한 메시가 아닙니다. 48 × 48에서 0.11%로 훌륭해 보이고, 96 × 96에서 1.09%로 나빠 보입니다. 나빠지는 게 아닙니다. 요소 2,300개쯤에서 공표값을 지나쳐 계속 간 것입니다. 가장 조밀한 네 메시에 w(n) = w∞ − C n^(−p)를 맞추면 ShellDKGQ의 w∞는 공표 급수값과 0.08% 이내, ShellMITC4는 1.57% 떨어져 있습니다. 물리와 일관됩니다. MITC4는 횡전단을 포함하고 DKGQ는 포함하지 않으며, 기준값은 Kirchhoff 얇은 쉘 급수해입니다. 실무적 결론은 단순합니다. 메시 하나가 공표값과 일치하는 것을 수렴으로 받아들이지 마십시오. 한 번 더 조밀하게 하고 어느 쪽으로 움직이는지 보십시오.

3. 평면요소: 락킹과 그 반대

길이 10, 깊이 1, ν = 0, 선단 전단력, 깊이 방향 1개 길이 방향 10개로 분할한 캔틸레버. ν = 0이면 이 단면에 대해 Timoshenko 해가 정확합니다.

w = P L³ / (3 E I) + 6 P L / (5 G A) = 0.040240

요소선단변위오차정체

tri310.00908−77.44%등변형률 삼각형

quad0.02680−33.40%쌍선형·완전적분 — 전단락킹

bbarQuad0.03440−14.51%b-bar, 체적 성분만

enhancedQuad0.04010−0.35%추가 가정변형률

SSPquad0.08000+98.81%적분점 1개, 안정화

전단락킹이 순수 quad의 −33%입니다. 쌍선형 요소는 휨 변형장의 곡률을 표현하지 못합니다. 그 변위장이 가짜 전단변형률을 강제하고, 그 가짜 전단이 휨으로 갔어야 할 에너지를 흡수합니다. 요소가 조금 뻣뻣한 게 아니라 3분의 1만큼 뻣뻣하고, 뻣뻣한 쪽에서 수렴합니다.

SSPquad의 +99%는 반대 방향 실패입니다. 적분점 하나로는 휨 모드를 아예 볼 수 없어 인위적으로 안정화하는데, 세장비 10에서는 그 안정화가 부족해 요소가 두 배 물러집니다. 두 실패 모두 경고 없이 깔끔하게 수렴합니다.

enhancedQuad는 처음부터 맞습니다. 요소 10개에 0.35%. 추가 가정변형률 정식화는 휨을 표현하려고 내부 자유도를 더한 것이고, 모델 수준에서는 비용이 없습니다. 첫 인자만 다른 같은 명령입니다.

요소1040160640

tri31−77.44%−46.51%−17.91%−5.17% (삼각형 1,280)

quad−33.40%−11.20%−3.06%−0.78%

bbarQuad−14.51%−4.14%−1.07%−0.27%

SSPquad+98.81%+14.03%+3.17%+0.77%

enhancedQuad−0.35%−0.16%−0.04%−0.01%

quad는 0.78%에 도달하려고 요소 640개를 쓰고, enhancedQuad는 10개로 그것을 이깁니다. 요소 64배를 치르고 더 나쁜 답. 그리고 tri31은 삼각형 1,280개로도 여전히 5% 벗어나 있습니다. 등변형률 삼각형은 휨 문제에 설 자리가 없고, 변형률 구배가 작은 까다로운 형상을 메우는 데에만 씁니다.

4. 솔리드: 같은 선택, 5천 배의 비용

똑같은 캔틸레버를 8절점 솔리드로. ν = 0이면 평면응력과 3차원 탄성이 일치하므로 목표는 같은 0.040240이고, 순수 stdBrick이 평면응력 quad 결과를 모든 자리까지 재현합니다. 두 메시가 같은 일을 하고 있다는 유용한 내부 검사입니다.

요소10806405,1205,120에서의 시간

stdBrick−33.40%−11.20%−3.06%−0.78%4,288 ms

bbarBrick−7.75%−3.56%−1.01%−0.26%4,468 ms

SSPbrick−0.35%−0.16%−0.04%−0.01%3,846 ms

요소 10개의 SSPbrick은 0.8 ms가 걸리고, 4.3초가 걸리는 요소 5,120개의 stdBrick보다 두 배 정확합니다. SSPbrick이 SSPquad와 전혀 다르게 거동한다는 점에 주의하십시오. 이 세장비에서는 안정화가 충분했고, 수치는 enhancedQuad와 정확히 같습니다. 교훈은 "SSP를 써라"도 "enhanced를 써라"도 아닙니다. 연속체 모델이라면 어느 쪽을 믿기 전에 같은 문제를 요소 두 가지로 돌려보라는 것입니다. 그 불일치는 재는 데 비용이 들지 않고, 여러분의 오차를 가둡니다.

5. 지반: 정확히 검사할 수 있는 주상체

30 m 지반 주상체, Vs = 200 m/s, ρ = 2.0 Mg/m³이므로 G = ρVs² = 80,000 kPa. 평면변형 요소를 한 열로 쌓고, 바닥을 고정하고, 각 높이의 양쪽 절점을 equalDOF로 묶어 순수 전단으로 변형하게 합니다. 균질 전단보의 주기는 정확해가 있습니다.

T_n = 4H / ((2n − 1) Vs) -> 0.6000, 0.2000, 0.1200, 0.0857 s

고유치에 들어가기 전에, 틀리면 ρ배를 치르는 측정 하나. SSPquad의 체적력 인자는 가속도가 아니라 단위체적당 힘입니다.

b2에 넣은 값바닥 요소의 연직응력정확해

-9.80665 (중력가속도)−289.296 kPa−578.592 kPa

-RHO*9.80665−578.592 kPa

단위중량 대신 중력가속도를 넣으면 연직응력이 정확히 ρ배 작게 나오고, 경고는 없으며, 그 뒤로 구속압에 의존하는 모든 지반 물성이 틀립니다.

모드주기종류닫힌해오차

10.600010전단0.6000000.002%

20.320719압축0.3207130.002%

30.200029전단0.2000000.014%

40.120048전단0.1200000.040%

50.106920압축0.1069040.014%

60.085782전단0.0857140.079%

요소 80개에서 모든 주기가 이론과 0.08% 이내이므로 모델은 맞습니다. 하지만 2차 모드는 압축 모드로 4H/Vp = 0.3207 s이지, 0.2000 s인 2차 전단모드가 아닙니다. 모드 번호로 "2차 모드"를 읽으면 주기가 60% 길게 나옵니다. 8부의 경고가 새 무대에 다시 나온 것입니다. 모드는 번호가 아니라 모드형상으로 분류하십시오. 그렇게 하는 두 줄입니다.

v = ops.nodeEigenvector(node, mode) # [ux, uy] horizontal = v[0]**2 / (v[0]**2 + v[1]**2) # 1.0000 전단, 0.0000 압축

6. 비선형 지반: 입력이 아니라 강도가 정한다

ElasticIsotropic을 점착력 40 kPa인 PressureIndependMultiYield로 바꾸고, 바닥을 Tabas FN 기록으로 0.045 g에서 0.9 g까지 스케일링해 흔듭니다. 탄성 주상체가 대조군입니다.

입력 PGA탄성c = 40 kPa

지표 PGA최대 전단응력지표 PGA최대 전단응력최대 전단변형률

0.045 g0.1458 g42.8 kPa0.0779 g23.7 kPa0.06%

0.180 g0.5834 g171.1 kPa0.1283 g43.9 kPa1.20%

0.360 g1.1667 g342.2 kPa0.1790 g45.4 kPa8.66%

0.900 g2.9168 g855.5 kPa0.2377 g46.2 kPa40.3%

입력이 20배 늘 때 지표에서는 3.1배 늘어납니다. 전단응력은 46.2 kPa에서 멈춥니다. 점착력의 1.15배이고, 그 약간의 초과는 항복면의 팔면체 정의와 감쇠응력 때문이며, 입력이 다섯 배로 커지는 동안 그 자리에 머뭅니다. 구성식이 약속한 그대로를 하는 것이고, 지반 모델에서 쓸 수 있는 가장 강한 검사입니다. 동원된 전단응력이 여러분이 준 강도를 넘어서는 안 됩니다. 한편 탄성 주상체는 855 kPa을 보고합니다. 자기가 들어본 적도 없는 강도의 21배이고, 0.9 g 입력에서 증폭 3.24배입니다. 그리고 모든 스텝에서 완벽하게 수렴합니다.

7. 시간 간격은 파동이 정하고, 회복 개입 횟수가 틀렸다고 알려준다

위 표의 첫 버전은 틀렸고, 모델은 제가 하마터면 무시할 뻔한 방식으로 그렇다고 말하고 있었습니다. 비선형 주상체를 Δt = 0.005 s로 — 7부와 9부의 골조 모델에서는 아무 문제 없던 간격으로 — 돌렸을 때, 해석은 9부 캐스케이드의 수렴 회복 개입을 143번 필요로 했고 그 전부가 성공했습니다.

Δt지표 PGA최대 전단응력최대 전단변형률회복 개입 횟수

0.005 s4.103 g43.83 kPa2.69%143

0.002 s0.172 g45.05 kPa7.13%0

0.001 s0.168 g45.05 kPa7.14%0

0.0005 s0.167 g45.05 kPa7.14%0

성긴 간격의 해석은 지표 가속도를 24배 과대평가하고 최대 전단변형률을 2.7배 과소평가했으며, 끝까지 완주했습니다. 9부에서 회복 개입 횟수는 결과 옆에 적어야 한다고 주장했는데, 여기가 그것을 증명하는 사례입니다. 143번의 개입은 간격이 너무 크다고 모델이 말하고 있었던 것이고, 그럼에도 나온 답은 24배 틀렸습니다.

이유는 수치가 아니라 물리입니다. 골조 모델의 시간 간격은 관심 있는 주기가 정합니다. 파동전파 모델의 시간 간격은 전단파가 요소 하나를 통과하는 시간이 정합니다. Vs = 200 m/s에 요소 1 m면 5 ms — 실패한 바로 그 간격입니다. 출발점으로 쓸 만한 규칙은 Δt ≤ h / (5 Vs)이고, 여기서는 1 ms이며, 그다음 간격을 반으로 줄여 답이 움직이지 않는지 확인하십시오.

8. 연속체 모델에서 확인할 것

ndf가 요소에 맞는가 — 평면 2, 솔리드 3, 쉘 6?

재료가 nDMaterial이고, 그 타입 문자열(PlaneStress, PlaneStrain)이 의도한 것인가?

체적력이 가속도가 아니라 단위체적당 힘인가?

같은 문제를 두 번째 요소 타입으로도 돌렸는가, 그리고 둘이 얼마나 벌어지는가?

메시를 한 번 더 조밀하게 했는가, 그리고 답이 움직였는가?

모델 어딘가에 닫힌해 검사가 있는가 — 주기, 정적 처짐, 알려진 깊이의 응력?

파동 문제라면 Δt가 요소 크기에 비해 충분히 작은가, 그리고 반으로 줄이면 답이 바뀌는가?

비선형 재료라면 동원된 응력이 지정한 강도 안에 있는가?

그 해석은 회복 개입을 몇 번 필요로 했고, 그 숫자가 보고되었는가?

9. 같은 모델을 다섯 단계 깊이로

입문자는 메시를 만들고 숫자를 읽습니다. 주니어는 숫자가 움직이지 않을 때까지 메시를 조밀하게 합니다. 실무 해석자는 quad가 휨에서 락킹한다는 것을 알고, 메시를 늘리기 전에 enhancedQuad나 SSPbrick을 꺼냅니다. 시니어는 정식화 두 가지를 돌려 그 불일치를 오차막대로 취급하고, 시간 간격을 습관이 아니라 요소 크기와 파속에서 정합니다. 전문가는 동원된 응력을 구성식 강도와 대조하고 모드를 형상으로 분류합니다. 틀린 모델이 흉내 낼 수 없는 두 가지가 그것이기 때문입니다.

다음 글은 이 모든 것을 보는 이야기입니다. 후처리와 시각화 — 변형형상, 모드형상, 파이버 응력 분포, 애니메이션 — 그리고 그림과 숫자가 어긋날 때 어떻게 할 것인가.

전체 스크립트

같은 10요소 메시 위의 평면요소 다섯 가지를 Timoshenko 해와 대조합니다.

# Part 11 - shear locking: five plane elements on the same ten-element mesh # units: kN, m import openseespy.opensees as ops L, H, TH, E, NU, P = 10.0, 1.0, 1.0, 1.0e5, 0.0, 1.0 I, A = TH * H**3 / 12.0, TH * H Gm = E / (2.0 * (1.0 + NU)) W_TIMO = P * L**3 / (3.0 * E * I) + 6.0 * P * L / (5.0 * Gm * A) def tip(ele_type, nx, ny): ops.wipe(); ops.model('basic', '-ndm', 2, '-ndf', 2) ops.nDMaterial('ElasticIsotropic', 1, E, NU) idx = {}; tag = 1 for i in range(nx + 1): for j in range(ny + 1): ops.node(tag, i * L / nx, j * H / ny - H / 2); idx[(i, j)] = tag; tag += 1 for j in range(ny + 1): ops.fix(idx[(0, j)], 1, 1) e = 1 for i in range(nx): for j in range(ny): n1, n2 = idx[(i, j)], idx[(i + 1, j)] n3, n4 = idx[(i + 1, j + 1)], idx[(i, j + 1)] if ele_type == 'tri31': ops.element('tri31', e, n1, n2, n3, TH, 'PlaneStress', 1); e += 1 ops.element('tri31', e, n1, n3, n4, TH, 'PlaneStress', 1); e += 1 elif ele_type == 'bbarQuad': ops.element('bbarQuad', e, n1, n2, n3, n4, TH, 1); e += 1 elif ele_type == 'SSPquad': ops.element('SSPquad', e, n1, n2, n3, n4, 1, 'PlaneStress', TH); e += 1 else: ops.element(ele_type, e, n1, n2, n3, n4, TH, 'PlaneStress', 1); e += 1 ops.timeSeries('Linear', 1); ops.pattern('Plain', 1, 1) for j in range(ny + 1): ops.load(idx[(nx, j)], 0.0, -P / ny * (0.5 if j in (0, ny) else 1.0)) ops.constraints('Plain'); ops.numberer('RCM'); ops.system('UmfPack') ops.test('NormDispIncr', 1e-12, 20); ops.algorithm('Linear') ops.integrator('LoadControl', 1.0); ops.analysis('Static'); ops.analyze(1) return -sum(ops.nodeDisp(idx[(nx, j)], 2) for j in range(ny + 1)) / (ny + 1) print(f'Timoshenko tip deflection {W_TIMO:.6f} m\n') print(f"{'element':14s}{'w (m)':>11s}{'error':>10s}") res = {} for et in ('tri31', 'quad', 'bbarQuad', 'enhancedQuad', 'SSPquad'): w = tip(et, 10, 1); res[et] = w print(f'{et:14s}{w:11.5f}{100*(w/W_TIMO - 1):9.2f}%') assert abs(res['enhancedQuad'] / W_TIMO - 1) < 0.01, 'enhancedQuad should be within 1%' assert res['quad'] / W_TIMO < 0.75, 'quad should lock badly on this mesh' print('\nOK the same ten elements span -77% to +99%; only the formulation changed')

실행하면 이렇게 나옵니다.

Timoshenko tip deflection 0.040240 m element w (m) error tri31 0.00908 -77.44% quad 0.02680 -33.40% bbarQuad 0.03440 -14.51% enhancedQuad 0.04010 -0.35% SSPquad 0.08000 98.81% OK the same ten elements span -77% to +99%; only the formulation changed

OpenSees 시리즈 전체

각 편은 닫힌해·독립적인 풀이·OpenSees 소스 중 하나와 대조해 검증했다. 순서대로 읽도록 썼지만 각 편이 자기 전제를 스스로 밝힌다.

1부 — 도메인 모델과 해석 조립

2부 — 단면·기하변환·분포하중

3부 — 파이버 단면·모멘트–곡률·푸시오버

4부 — 질량·감쇠·시간적분

5부 — 라이브러리 읽기와 재현 가능한 워크플로

6부 — 철근콘크리트 파이버 단면과 구속

7부 — 실제 지진기록·층간변위·증분동적해석

8부 — 3차원: 강막·비틀림·leaning column

9부 — 수렴하지 않을 때: 실패·알고리즘·한계점

10부 — Tcl 읽기와 공개 스크립트 이식

11부 — 쉘·솔리드·지반: 락킹·메시·부지응답 — 지금 보는 글

12부 — 후처리: recorder·질의, 그리고 그림이 증명하는 것

13부 — 좌굴과 한계점: 분기점·스냅백·호장법

14부 — Windows: 설치·파일 형식·조용한 실패

15부 — 검증된 정답이 딸린 연습문제와 한영 용어집

16부 — 어떻게 동작하는가: 시행상태·확정상태·해석 루프

17부 — 커뮤니티가 말하는 것, 측정해보면

18부 — 절점 변위에서 파이버 응력까지

19부 — 면진과 감쇠: 속도 의존 요소와 베어링

참고 자료

OpenSees, PinchedCylinder.tcl — 2절의 벤치마크이자 기준값의 출처.

Lindberg, G. M., Olson, M. D., Cowper, G. R., "New Developments in the Finite Element Analysis of Shells", Quarterly Bulletin of the Division of Mechanical Engineering and the National Aeronautical Establishment, NRC Canada, vol. 4, 1969.

OpenSeesPy, quad, bbarQuad, enhancedQuad, SSPquad, tri31.

OpenSeesPy, stdBrick, bbarBrick, SSPbrick.

OpenSeesPy, ShellMITC4, ShellDKGQ, PlateFiber 단면.

OpenSeesPy, PressureIndependMultiYield, updateMaterialStage — 6절의 2단계 지반 워크플로.

Kramer, S. L., Geotechnical Earthquake Engineering, Prentice Hall, 1996 — 균질 전단보와 1차원 부지응답.

Elgamal, A., Yang, Z., Parra, E., "Computational modeling of cyclic mobility and post-liquefaction site response", Soil Dynamics and Earthquake Engineering 22(4), 2002 — 다중항복면 재료.

출처 확인 2026-08-31. 모든 결과는 OpenSeesPy 3.8.0에서 산출했으며 본문의 스크립트로 재현할 수 있습니다. 요소 순위는 두 개의 특정 문제 — 휨을 받는 세장한 캔틸레버와 얇은 쉘 — 에 대한 측정치이며 다른 변형 모드가 지배하는 문제로 옮겨지지 않습니다. 옮겨 쓸 수 있는 결론은 두 정식화의 불일치가 문제가 될 만큼 크다는 것이고, 그래서 둘 다 돌려야 한다는 것입니다. 시간 측정은 한 대의 기계에서 나온 것으로 비율로 읽어야 합니다. 결과는 프레임워크의 거동을 보여주는 것으로, 프로젝트별 해석이나 설계기준 검토를 대체하지 않습니다.