Advanced
Development of Particle Simulation Method for Analysis of Fluid-Structure Interaction Problems
Development of Particle Simulation Method for Analysis of Fluid-Structure Interaction Problems
Journal of Ocean Engineering and Technology. 2013. Apr, 27(2): 53-58
Copyright © 2013, Korean Society of Ocean Engineers
  • Received : January 29, 2013
  • Accepted : April 19, 2013
  • Published : April 30, 2013
Download
PDF
e-PUB
PubReader
PPT
Export by style
Article
Author
Metrics
Cited by
TagCloud
About the Authors
성철 황
종천 박
창용 송
영훈 김

Abstract
Recently, some fluid-structure interaction (FSI) problems involving the fluid impact loads interacting with structures, such as sloshing, slamming, green-water, etc., have been considered, especially in the ocean engineering field. The governing equations for both an elastic solid model and flow model were originally derived from similar continuum mechanics principles. In this study, an elastic model based on a particle method, the MPS method, was developed for simulating the FSI problems. The developed model was first applied to a simple cantilever deflection problem for verification. Then, the model was coupled with the fluid flow model, the PNU (Pusan National University modified)-MPS method, and applied to the numerical investigation of the coupling effects between a cantilever and a mass of water, which has variable density, free-falling to the end of the cantilever.
Keywords
1. 서 론
유체-구조 연성(Fluid-structure interaction, FSI)이란 수치적으로 유체의 유동해석과 고체의 구조해석에 있어 상호작용을 고려하는 것이다. 최근 항공( Park et al., 2011 ), 조선·해양( Koo et al., 2009 ; Lee and Rlee, 2011 ; Shin, 2007 ,: Shin, 2011 ) 및 의료( Lee et al., 2011b ) 등 다양한 분야에서 FSI 해석에 관한 연구가 진행 중이다. 특히, 조선·해양 분야에서는 선박 또는 해양 구조물이 운행 중 갑판침수(Green water), 슬로싱(Sloshing), 또는 슬래밍(Slamming)과 같은 원인에 의해 유체 충격력을 받게 되면, 선체 구조는 변형이나 응력이 발생하게 되는데, 이때 발생한 응력과 변형은 이들의 안전성과 운항성에 영향을 미치게 된다. 그러나 이전의 수치 해석 기법은 유동해석과 구조해석을 각각 독립적으로 수행해 왔으며 상호작용에 대한 연성을 포함하지 않는 경우가 대부분이었다. 그러나 점차 선박의 대형화 및 해양구조물의 개발이 가속화됨에 따라 정도 높은 해석과 유체-구조 상호연성 해석의 필요성이 강조되고 있다.
유체-구조 상호작용 연성해석이란 유동해석과 구조해석의 연동해석을 의미한다. 연동 해석은 먼저 유동 해석을 통해 구조물의 표면에 작용하는 유체 충격력을 계산하고, 구조 해석시 외력에 대한 하중조건으로 부여하여 구조물의 변형과 응력을 계산한다. 계산된 구조물의 변형 및 변형속도는 다음의 유동 해석에서 역시 경계조건이 되어 참조가 되는 일련의 과정을 포함하게 된다. 일반적으로, 탄성체 모델의 지배방정식은 선형 탄성체의 평형 방정식을, 유동 해석에는 Navier-stokes(N-S) 방정식과 연속 방정식을 사용한다. 이러한 해석기법은 최근 상용 해석 프로그램들을 중심으로 다양한 분야에서 연동해석의 지원을 시작하고 있다( Choi et al., 2009 ; Kennedy and Solbrekken, 2011 ).
유동 해석의 경우, 지금까지 연산 속도와 계산기 용량의 한계 등으로 효율적인 수치 연산에 유리한 오일러(Eulerian) 관점의 격자기반 해석 방법들이 주로 개발되어 왔다. 하지만, 최근 컴퓨터의 급속한 발전으로 대용량의 수치 연산이 가능해져 이로 인해 기존의 수치연산에서 배제되어 오던 라그란지(Lagrangian) 관점의 해석 기법인 이동 입자 기반의 시뮬레이션 기법들이 제안되고 있다.
Monaghan(1988) 이 제안한 SPH(Smooth particle hydro-dynamics)법과 Koshizuka and Oka(1996) 에 의해 제안된 MPS(Moving particle simulation)법은 대표적인 입자를 이용한 유동해석 기법이다. 격자법에 비해 입자법은 형상의 재현이 쉽고 격자 생성이 손쉬우며 특히 물질의 분해 및 결합 등과 같은 극도로 비선형성이 강한 물리적 현상을 표현하는데 큰 장점을 갖고 있다. 그러나 입자법의 경우 비물리적인 압력 진동과 막대한 계산 속도는 해결해야 할 단점으로 지적되고 있다.
본 연구에서는 입자법인 MPS법( Kondo et al., 2005 )을 이용하여 탄성체 해석 모델을 개발하고 유체와의 연성 해석에 대한 가능성을 확인하였다. 먼저, 탄성체 모델을 외팔보(Cantilever beam)의 처짐 문제에 적용하여 해석해와의 비교를 통해 개발된 모델의 타당성을 검증하였다. 또한, 유동해석에 MPS법( Koshizuka and Oka, 1996 )을 개량한 PNU-MPS법( Lee et al., 2011a )을 이용하여 탄성체 모델과의 연성 모델을 개발하고 외팔보의 끝단에서 자유낙하 하는 정사각형의 물덩이와의 유체-구조 상호작용 문제에 적용해 보았다.
2. 수치해석 기법
- 2.1 지배 방정식
선형 탄성체의 평형 방정식은 다음과 같다.
PPT Slide
Lager Image
여기서, ui 는 변위 벡터이고, ρ 는 밀도, t 는 시간, 𝜖 kk 는 체적변형률 텐서, 𝜖 ij 는 변형률 텐서, fi 는 외력, 𝛿 ij 는 크로넬커 델타(Kronecker's delta), 그리고 λ와 𝜇는 라메(Lame) 상수이며 다음과 같다.
PPT Slide
Lager Image
PPT Slide
Lager Image
단, E 는 탄성계수(Young’s modulus)이고, ν 는 포아송 비(Poisson ratio)이다.
- 2.2 이산화
선형 탄성체의 평형 방정식의 지배방정식에 관한 이산화에는 MPS법( Kondo et al., 2005 )의 연산모델을 사용하기로 한다.
식 (1)의 첫 번째 항은 다음과 같이 표현 할 수 있다.
PPT Slide
Lager Image
여기서 주변 입자 j 와 입자 i 의 체적변형률( ϵkkδij ) ij 는 다음과 같이 정의한다.
PPT Slide
Lager Image
단, rij 는 입자 i j 의 상대 위치 벡터이며,
PPT Slide
Lager Image
는 초기의 상대 위치 벡터이다. 이를 이용하여 ( ϵkk ) i 는 다음과 같이 표현된다.
PPT Slide
Lager Image
여기서, d 는 계산 공간의 차원, Wij 는 커넬함수이며 본 연구에서는 Koshizuka and Oka(1996) 에 의해 제안된 다음 식을 사용하였다.
PPT Slide
Lager Image
단, re 는 입자간 상호작용 반경이다. 식 (6)의 n0 는 초기 입자수 밀도로 다음과 같이 초기 주변 입자의 커넬함수의 합이다.
PPT Slide
Lager Image
식 (6)으로부터 구한 ( ϵkkδij ) i 를 이용하여 식 (4)는 MPS 법의 발산모델을 적용하여 다음과 같이 이산화 할 수 있다.
PPT Slide
Lager Image
여기서, W 0 는 초기 상대위치에 대한 커널함수의 값이다.
다음으로 식 (1)의 우변 두 번째 항은 전단변형률에 관한 항으로 다음과 같이 정의 할 수 있다.
PPT Slide
Lager Image
식 (10)에서 전단 변형률은 변위벡터의 라플라스(Laplacian) 연산 결과와 같아지며, MPS법의 라플라스 모델을 이용하면 다음과 같이 이산화 할 수 있다.
PPT Slide
Lager Image
식 (11)의 변위 벡터 차는 다음과 같이 상대 위치 벡터의 차로서 구할 수 있다.
PPT Slide
Lager Image
식 (11)와 같이 상대 위치 벡터에 의한 전단변형률 성분은 초기 상대 위치로 되돌아가려는 성질을 가지게 된다. 그러나 탄성체가 강체 회전을 할 경우 강체 회전에 의한 상대위치 변화성분은 제외되어야 한다. 따라서 강체회전에 의한 성분이 제외된 수정된 상대위치 벡터는 아래와 같다.
PPT Slide
Lager Image
여기서 Ri 는 입자 i 의 회전에 대한 회전 연산자이며 초기 상대 위치 벡터를 입자 i 와 입자 j 의 상대위치에 회전 변환시킨 평균치로서 정의한다. 식 (13)을 식 (11)에 대입하여 최종 전단변형률 성분을 다음과 같이 얻을 수 있다.
PPT Slide
Lager Image
최종 변위 벡터의 가속도 성분은 식 (9)와 식 (14)의 합으로 구성되지만, 여기에는 탄성체의 강체 회전성분이 고려되지 않았다. 따라서 추가적으로 각 운동량 보존법칙을 이용하여 각 입자들의 강체회전성분에 대한 연산을 수행할 필요가 있다. 이때 각운동량 보존식은 다음과 같다.
PPT Slide
Lager Image
여기서 ω 는 각속도, r 은 회전 중심으로부터의 거리벡터, F 는 전단력이며 전단변형률에 의한 가속도항으로 대체할 수 있다. 또한 I 는 입자의 질량 관성 모멘트이고 다음과 같이 주어진다.
PPT Slide
Lager Image
여기서 m 은 입자의 질량, l0 는 초기 입자간 거리이며 각 입자의 가로세로의 길이를 의미한다.
따라서 식 (15)는 입자 i 에 대해 다음과 같이 이산화 된다.
PPT Slide
Lager Image
이렇게 구해진 각속도 변화량과 속도 변화량을 이용하여 아래의 식들과 같이 다음 시간 단계의 위치와 속도, 각속도, 그리고 회전연산자를 구하게 된다.
PPT Slide
Lager Image
PPT Slide
Lager Image
PPT Slide
Lager Image
PPT Slide
Lager Image
여기서,
PPT Slide
Lager Image
은 회전연산을 나타낸다.
- 2.3 연동(Coupling) 기법
비압축성 유동장의 지배방정식은 다음과 같은 연속방정식과 N-S 방정식이다.
PPT Slide
Lager Image
PPT Slide
Lager Image
여기서 u 는 속도이고, 𝜇는 점성계수, p 는 압력, ρ 는 밀도, g 는 중력가속도이다.
식 (22)와 (23)은 PNU-MPS법( Lee et al., 2011a )의 입자간 상호작용 모델을 이용하여 시뮬레이션을 수행한다. 유동해석의 이산화방법은 참고문헌( Lee et al., 2011a )에 상세히 나타나 있다.
Fig. 1 에는 유체-구조 상호작용에 관한 개략적인 연성해석법을 나타내며, 다음과 같이 정리할 수 있다.
PPT Slide
Lager Image
Diagram of FSI coupling
  • ① PNU-MPS을 이용하여 유동해석을 수행하여 각 벽면상에 작용하는 압력을 계산(이때 유동해석의 벽면경계입자들은 강체로 적용하며, 구조해석으로부터 구한 구조입자의 변형과 변형속도를 이용하여 경계조건으로 부여)
  • ② 유동해석을 통해 구한 각 벽면경계입자들에 작용하는 압력에 각 입자의 미소면적을 곱하여 식 (1)의 우변 두 번째 항인fi=PiΔAi를 산출하여 구조해석의 외력에 대한 하중조건으로 적용(단, 이때 각 입자의 미소면적은
  • ③ 구조해석을 수행하여 동일한 시각에서 각 구조입자들의 변형 및 변형속도 계산
이로써 한 시간 단계의 계산이 완료되고 일정한 시간까지 반복하여 연동해석을 수행한다. 이때 시간적분 방법은 시차제법(Staggered method)을 사용하여 매 시간 양방향 연동해석을 수행하였다.
3. 시뮬레이션 결과
- 3.1 외팔보의 처짐 문제
2.1 절에서 설명한 탄성체 해석 기법을 이용하여 2차원 외팔보의 처짐 문제에 적용하였다. Fig. 2 는 시뮬레이션의 초기 조건을 나타난다. 외팔보의 길이 L =0.16m, 두께 T =0.02m인 외팔보의 오른쪽 끝단에서는 외력 P =0.005 N 이 작용하며, 왼쪽 고정단은 벽면에 완전히 고정된 것으로 가정한다. 여기서 E 는 1.0×10 5 Pa이고, ν 는 0.3이다. 입자간 간격은 0.005m이며 시뮬레이션에 사용된 입자는 총 180개이다. 가중치 함수의 임계거리 re 는 초기 입자배치 간격의 2.1배를 사용하였다.
PPT Slide
Lager Image
Initial setup for cantilever deflection problem
Fig. 3 는 시뮬레이션에서 구한 외팔보의 처짐량을 식 (24)의 해석해( Ugural and Fenster, 2011 )와 비교한 결과이다. 이 때 해석해는 세장비가 작은 외팔보 문제에 대한 해석해로 포아송비에 의한 전단 변형이 고려된 식이다. 전체적으로 보의 처짐에 대한 경향이 잘 일치하는 것을 알 수 있다. 이 때, 최대 처짐에 대한 상대 오차는 3.3% 이내이다. 상대 오차의 발생은 식 (24)가 연속체역학에 기반하여 산출된 결과인 반면, 본 연구에서 적용한 시뮬레이션 기법과 같은 이산화 수치해석법 상의 입자구성 과정 중에 발생하는 상대적 수치 오차에 기인한 것으로 볼수 있다. Fig. 4 는 외팔보에 작용하는 응력 분포를 나타내며 식(25)의 해석해와 비교하였다. 여기서 Z 는 단면 계수이다. 시뮬레이션 결과는 해석해와 비교하여 고정단에서 응력의 최대치를 보이며 외팔보의 끝단으로 가면서 0에 수렴하는 경향이 잘 일치하는 반면, 응력의 크기에는 다소 차이를 보인다. 이 때, 각 지점에서 해석해에 대한 평균 오차는 약 8%이다. 이러한 오차는 입자구성 상의 이산화에 따른 수치적 오차에 기인하며 향후 최적화된 입자구성 모델의 정도 향상 등을 통해 개선이 가능할 것으로 보인다. 본 시뮬레이션은 Intel(R) Core(TM) i7 CPU 2.8GHz, 4GHz Memory의 PC에서 수행되었으며 총 CPU 시간은 965초가 소요되었다.
PPT Slide
Lager Image
Comparison of deflection of cantilever
PPT Slide
Lager Image
Comparison of Stress distribution of cantilever
PPT Slide
Lager Image
PPT Slide
Lager Image
- 3.2 FSI 연동 해석
본 절에서는 탄성체 해석 기법을 유동해석 기법인 PNU-MPS법과 연동하여 FSI 해석에 관한 시뮬레이션을 수행하였다. 시뮬레이션의 초기조건은 Fig. 5 와 같다. 단, 외팔보의 조건은 3.1절과 동일하다. 외팔보의 끝단에는 h =0.1(m) 높이에서 xwater × ywater =0.02×0.02(m 2 )의 정사각형 물덩이가 자유낙하 하여 구조물과 상호작용하는 것으로 설정하였다. 구조물에 작용하는 힘을 다양하게 시험하기 위해 유체의 밀도 ρ 를 500, 1000, 1500(kg/m 3 )으로 변화시켜 보았다. 이때 초기 입자간 간격은 0.002m이고 사용된 입자의 수는 외팔보가 830개, 유체가 100개이다. 한편, 탄성체 해석의 경우 시간 간격 Δt 는 10 −6 (sec)이고 유동해석의 경우 Courant 수를 0.2로 유지하도록 가변적인 시간간격을 사용하였다. 따라서 본 시뮬레이션의 경우 Δt =10 −3 ~2×10 −4 (sec)의 범위가 사용되었다.
PPT Slide
Lager Image
Initial setup for FSI coupling simulation
Fig. 6 ρ =1000(kg/m 3 )일 때 자유낙하 하는 물덩이가 외팔보에 충돌하여 상호작용하는 모습의 시간 변화를 나타낸다. 이때 유체입자의 색은 각 입자의 압력을 나타내고 구조입자의 변형속도는 화살표로 표시하였다. 유체의 경우 공기의 저항 및 표면장력을 무시했으므로 초기조건의 형상 그대로 낙하하여 외팔보 끝단에 충돌하여 탄성체에 큰 충격력을 주게 된다. 유체에 의한 충격력에 의해 외팔보는 아래로 처진 후 다시 위로 올라오는 모습을 확인 할 수 있다. 이때, 탄성체에는 끝단으로부터 고정단 쪽으로 점차 동적외력이 전달되는 모습이 보인다.
PPT Slide
Lager Image
Distribution of particles on FSI problem.
Fig. 7 은 유체의 밀도 차에 의한 외팔보 끝단의 처짐 변위에 대한 시간 변화를 나타낸다. 밀도가 커짐에 따라 전체적인 경향은 다소 유사하지만 최대 처짐량이 증가하며 진동하는 것을 알 수 있다. 이때의 최대 처짐량을 Table 1 에 정리하여 비교하였다.
PPT Slide
Lager Image
Time-history of deflection at the end of cantilever with various water density
Comparison of maximum deflection the end of cantilever with various water density
PPT Slide
Lager Image
Comparison of maximum deflection the end of cantilever with various water density
본 시뮬레이션은 Intel(R) Core(TM) i7 CPU 2.8GHz, 4GHz Memory의 PC에서 수행되었고 0.5초까지의 총 CPU 시간은 2,527초가 소요 되었다.
5. 결 론
본 연구에서는 입자법을 이용하여 유체-구조 상호연성 해석을 위한 탄성체 해석 모델을 개발하였다. 개발된 탄성체 모델은 MPS법의 입자간 상호작용 모델( Kondo et al., 2005 )을 이용하여 이산화 되었다. 먼저, 모델의 타당성을 검증하기 위하여 단순한 외팔보 끝단에 집중하중이 작용할 경우의 처짐 문제에 적용하였다. 시뮬레이션 결과는 해석해와 정성적·정량적으로 잘 일치함을 알 수 있었다. 다음으로, 유동해석에 PNU-MPS법( Lee et al., 2011a )을 이용하여 탄성체 모델과의 연성 모델을 개발하고 다양한 밀도의 유체 덩어리가 자유낙하 하여 외팔보에 충돌하는 문제에 적용하여 보았다. 이를 통해, 입자법을 이용한 유체-구조 상호작용 해석 기법에 대한 가능성을 확인할 수 있었다. 향후 포아송 비를 포함한 3차원 처짐량에 대한 고려 및 이산화 과정의 고도화를 통해 다양한 유체-구조 해석 문제에 적용하여 개발된 탄성체 모델과 유체해석 기법간의 상호작용 연성 기법의 타당성 검증을 수행할 예정이다.
본 논문은 2011년도 정부(교육과학기술부)의 재원으로 한국연구재단의 지원을 받아 수행된 기초연구사업(2011-0011220)과 2012년도 정부(교육과학기술부)의 재원으로 한국연구재단의 지원을 받아 수행된 기초연구사업(2012R1A1A1002897)임.
References
Choi C.R. , Kim C.N. , Hong T.H. 2009 Blood Cell Dynamics in a Simple Shear Flow using an Implicit Fluid-Structure Interaction Method Based on the ALE Approach World Academy of Science, Engineering and Technology 49 685 - 690
Kennedy J.C. , Solbrekken G.L. 2011 Coupled Fluid Structure Interaction(FSI) Modelling of Parallel Plate Assembles Proc. ASME 2011 IMECE 2011 1 - 9
Kondo M. , Koshizuka S. , Suzuki Y. 2005 Application of Sysmplectic Scheme to Three-dimensional Elastic Analysis Using Mps Method Transactions of the Japan Society of Mechanical Engineers 72 65 - 71
Koo J. , Jeong W. , Cho S. , Ahn S. , Kim D. 2009 Study on Adaptive FSI Technique for Hydroelastic Analysis of Liquified Fuel-Storage Tank Proceedings of the Korea Society for Noise and Vibration Engineering Conference 496 - 497
Koshizuka S. , Oka T. 1996 Moving-particle Semi-implicit Method for Fragmentation of Incompressible Fluid Nuclear Science and Engineering 123 421 - 434
Lee B.H. , Park J.C. , Kim M.H. , Hwang S.C. 2011 Step-by-step Improvement of MPS Method in Simulating Violent Free-surface Motions and Impact- loads Comput. Methods Appl. Mech. Engrg. 200 1113 - 1125    DOI : 10.1016/j.cma.2010.12.001
Lee H. , Rhee S.H. 2011 Fluid-structure Interaction Analysis of Two-dimensional Flow Around a Moving Cylinder Proceeding of Korean Society of Computational Fluid Engineering 68 - 74
Lee S.H. , Kang S. , Hur N. 2011 Development of Program for the Fluid-structure Interation Analysis on Hemodynamics of CerebralArtery Proceeding of Korean Society of Computational Fluid Engineering 355 - 358
Monaghan J.J. 1988 An Introduction to SPH Comput. Phys. Comn. 48 89 - 96    DOI : 10.1016/0010-4655(88)90026-4
Park K.H. , Min J.K. , Kim J.K. , Kang S.H. , Kim S.J. , Park S.H. 2011 A Numerical Study on the Effect of a Microfin with a Flexible Up-down Movement on Heat Transter using a Fluid-structure Interaction (FSI) Method Journal of the Korea Society for Precision Engineering 28 (8) 975 - 983
Shin S. 2007 Numerical Simulation on Fluid-structure Interaction of a Two-dimensional Orbiting Flexible Foil Journal of Computational Fluid Engineering 12 (2) 37 - 45
Shin S. 2011 Simulation of Fluid-Structure Interaction of a Towed Body using an Asysmmetric Tension Model Journal of Computational Fluid Engineering 16 (1) 7 - 13
Ugural A.C. , Fenster S.K. 2011 Advanced Mechanics of Materials and Applied Elasticity Fifth Edition Prentice Hall