Research Paper

Journal of the Computational Structural Engineering Institute of Korea. August 2020. 209-216
https://doi.org/10.7734/COSEIK.2020.33.4.209


ABSTRACT


MAIN

  • 1. 서 론

  • 2. 본 론

  •   2.1 이론적 배경

  •   2.2 해석 대상

  •   2.3 연구 진행

  •   2.4 연구 결과

  • 3. 결 론

1. 서 론

다양한 용도와 주행 특성에 따라 세분된 모터사이클 중에서 투어러라 불리는 부류는 장거리 여행을 즐기기 위해 생겨난 기종들이다. 안락하게 즐길 수 있는 편안한 자세와 넉넉한 배기량에서 오는 높은 출력은 필수 요소라 할 수 있다. 스포츠 투어러는 그러한 투어러에서 파생된 장르로, 투어러의 편안한 자세는 유지하면서도 스포츠 모델 못지않은 운동성을 지닌 장르이다.

이러한 투어러 모델들은 장거리 여행을 상정하는 그 특성상 짐을 수납할 공간들을 확장하기 위하여 Fig. 1과 같이 사이드 케이스와 탑케이스를 설치하는 경우가 많고 대부분의 제조사에서 공식 혹은 서드파티 액세서리의 형식으로 지원되고 있다. 그런데 많은 기종에서 탑케이스를 설치 후 Figs. 2~3과 같이 주행 시 후방프레임이 파손되는 사건이 보고되고 있다(Kim, 2018).

http://static.apub.kr/journalsite/sites/jcoseik/2020-033-04/N0040330401/images/Figure_jcoseik_33_04_01_F1.jpg
Fig. 1.

FJR 1300 with top case

http://static.apub.kr/journalsite/sites/jcoseik/2020-033-04/N0040330401/images/Figure_jcoseik_33_04_01_F2.jpg
Fig. 2.

Rear frame of FJR 1300

http://static.apub.kr/journalsite/sites/jcoseik/2020-033-04/N0040330401/images/Figure_jcoseik_33_04_01_F3.jpg
Fig. 3.

Broken rear frame

그러나 제조사들은 130km/h 이상 주행을 자제하라는 권고 외엔 아직 뚜렷한 해법을 내놓지 않은 상태이다.

후방프레임이 파손될 경우 탑케이스에 기대어 있던 동승자가 크게 다칠 우려가 있다. 동승자가 없더라도 파손된 프레임으로 인해 발생하는 탑케이스의 진동은 차량의 흔들림을 유발하여 사고로 이어질 가능성이 크다. 이처럼, 프레임이 파손되는 것은 운전자의 안전을 심각하게 위협할 수 있다. 그런 점에서 후방프레임은 정해진 규격의 탑케이스와 적재 제한 무게를 준수했을 경우 가동 범위 내에서 항상 안전한 주행성능을 확보할 수 있어야 한다. 그러므로 후방프레임의 구조적 거동에 대한 해석 및 이해가 필요하다.

선행 연구에서 Park(1998)은 모터사이클 차체의 진동특성을 분석하는 기법을 소개했으며 Lee 등(2016)은 프레임의 구조 강성 강화를 위한 설계변수해석을 수행했다. 기존의 연구들은 NISA II나 ABAQUS와 같은 상용 해석 솔루션들을 이용하여 해석이 이루어졌다.

이에 본 연구에서는 상용 해석프로그램이 아닌 co-rotational (CR) shell 해석기법을 적용하여 자체 개발된 코드를 통해 후방프레임 형상에 대해 해석을 수행하였다. 먼저 유사하지만 단순한 형상에 대해 상용 소프트웨어와 CR 기법으로 해석한 결과를 비교하여 CR기법의 신뢰성과 우위를 검증하였고, 이후 SOLIDWORKS를 활용해 모터사이클 후방프레임을 3차원 쉘 요소로 모델링한 뒤 응력해석과 공진 주파수 해석을 수행하였다. 이후 프레임의 파손을 방지하기 위해 설계 변경을 위한 매개변수를 설정하여 해석을 수행하였다. 본 연구에서 CR 쉘 요소의 해석은 EDISON 프레임워크를 활용하여 실행하였다. EDISON은 과학기술정보통신부 주도로 만들어진 계산과학 공학 솔루션으로, 전산 열유체, 나노 물리, 계산화학, 전산 설계 등 다양한 분야에서 각 연구실이 개발한 해석기법들을 한국과학기술정보연구원의 클라우드 컴퓨팅 기능을 활용하여 실행할 수 있는 서비스로, 계산 수행자의 PC 환경에 구애받지 않고 high-fidelity 계산을 수행한다는 장점이 있다. CR기법과 비교 검증하기 위해 사용된 상용 프로그램으로는 MSC NASTRAN/ PATRAN과 ANSYS가 사용되었다.

2. 본 론

이 장에서는 이론적 배경들을 소개한 후, 해석 대상물에 대한 설명, 그리고 해석 과정과 결과를 서술하고자 한다.

2.1 이론적 배경

이번 장에서는 이번 연구를 수행하면서 사용한 CR 기법, 그리고 CR기법 해석에서 국부 요소로 사용된 OPT-DKT 쉘 요소를 소개한다. 이후 CR기법과 OPT-DKT 쉘 요소가 사용된 해석과 상용 소프트웨어를 통한 해석과의 비교를 간단한 예제를 통해 수행한다.

2.1.1 CR 기법

본 논문에서 사용한 기법은 co-rotational(CR)기법으로, 구조물의 기하학적 비선형을 고려하는 유한요소 해석기법의 일종이다. CR기법은 Fig. 4와 같이 구조물의 전체 거동을 CR 좌표계 기준으로 강체 거동과 순수 구조 변형으로 분리한 뒤, 정립된 좌표계를 따라 전역 좌표계의 내력 벡터와 각 요소의 강성 행렬을 결정하는 기법이다(Cho and Shin, 2015a). CR기법은 요소에 독립적으로 적용 가능하여 기존의 선형 요소들을 비선형해석으로 확장 가능하다는 장점이 있고, 대 변위의 작은 변형을 갖는 문제에 적합하다. CR 쉘 요소의 강성 행렬과 내력 벡터는 식 (1)과 같이 정의된다.

$$p=\frac{\partial U}{\partial u},\;K=\frac{\partial^2U}{\partial u^2}=\frac{\partial p}{\partial u}=K_M+K_G$$ (1)
http://static.apub.kr/journalsite/sites/jcoseik/2020-033-04/N0040330401/images/Figure_jcoseik_33_04_01_F4.jpg
Fig. 4.

Coordinates in CR formulation

이때 강성 행렬 K는 재료 강성 행렬과 기하학적 강성 행렬의 결합으로 식 (2), (3)과 같이 된다.

$$K_M=T^TK_{M0}T$$ (2)
$$K_G=T^TK_GT+K_G^\nu$$ (3)

이때, KM0는 국부 요소의 재료 강성 행렬이며, 자세한 정식화 과정은 참고문헌(Cho and Shin, 2015b)에 나타나 있다.

2.1.2 OPT-DKT 쉘 요소

본 연구에서는 상기한 CR기법의 국부 요소로 optimal triangle membrane and discrete Kirchhoff triangle(OPT-DKT) 쉘 요소(Joo et al., 2017)를 도입하였다.

OPT-DKT 쉘 요소는 OPtimal Triangle membrane(OPT) 요소와 Discrete Kirchhoff Triangle(DKT) 요소가 결합한 형태로, 각각 membrane과 bending 요소를 담당하게 된다. Fig. 5에 묘사된 것처럼 각 요소는 3개의 절점을 가지며 각 절점에 대해 변형 자유도와 회전자유도를 3개씩 갖게 된다. OPT-DKT 쉘 요소의 강성 행렬은 본 논문 대상인 프레임의 알루미늄 합금 재질과 같은 등방성 재료의 경우 식 (4)와 같이 나타나게 된다.

$$K=\begin{bmatrix}K_m{(u_i,v_i,z_i)}_{i=1,2,3}&0\\0&K_b{(u_i,v_i,z_i)}_{i=1,2,3}\end{bmatrix}$$ (4)
http://static.apub.kr/journalsite/sites/jcoseik/2020-033-04/N0040330401/images/Figure_jcoseik_33_04_01_F5.jpg
Fig. 5.

OPT-DKT shell element

이러한 요소는 얇은 구조물에서 발생할 수 있는 잠금 현상(locking effect)을 예방할 수 있다. 본 연구에서 다루고 있는 후방프레임의 경우 경량화를 위해 내부가 비어 있는 얇은 형상을 지니기 때문에 이러한 쉘 요소를 사용하여 해석을 진행하였다.

2.1.3 CR기법 검증 예제

CR기법을 사용하여 해석한 결과와 상용 소프트웨어를 사용하여 해석한 결과를 비교하기 위해 사용한 예제는 Fig. 6과 같이 등방성 재질의 쉘로 이루어진 형상이다. 치수는 mm단위로 표시되어 있으며, 두께는 2mm이다. 경계 조건으로는 먼저 좌측 하단 부분을 고정단으로 구속하였다. 또한, 화살표로 표시된 부분에 하중을 설정하였다. 하중은 비교의 편의성을 위해 500N으로 설정하였다. 재료로는 본 형상과 동일하게 알루미늄 합금 6061-T6을 사용하였으며, 해당 합금의 물성치는 Table 1에 실려 있다.

http://static.apub.kr/journalsite/sites/jcoseik/2020-033-04/N0040330401/images/Figure_jcoseik_33_04_01_F6.jpg
Fig. 6.

Sample geometry

Table 1.

Properties of Al6061-T6

Density 2.70kg/m3
Young's Modulus 68.9GPa
Tensile Strength 276MPa
Ultimate Strength 310MPa
Poisson Ratio 0.33

요소의 크기를 조절하여 수를 조절하였으며, 각각 0.05m부터 0.005m까지 조절하며 수행하였다. Fig. 7은 이러한 과정을 거쳐 만든 이산화 형상이다. 수행 결과 정적 해석을 통해 얻은 변위는 Fig. 8과 같았다. Fig. 9에 나타난 MSC NASTRAN을 통해 계산한 변위 분포와 상당히 비슷한 결과를 보였다.

http://static.apub.kr/journalsite/sites/jcoseik/2020-033-04/N0040330401/images/Figure_jcoseik_33_04_01_F7.jpg
Fig. 7.

Discretized geometry of sample model obtained by CSD pre

http://static.apub.kr/journalsite/sites/jcoseik/2020-033-04/N0040330401/images/Figure_jcoseik_33_04_01_F8.jpg
Fig. 8.

Displacement by CR shell analysis(Mesh element size: 0.005m)

http://static.apub.kr/journalsite/sites/jcoseik/2020-033-04/N0040330401/images/Figure_jcoseik_33_04_01_F9.jpg
Fig. 9.

Displacement by NASTRAN analysis (Mesh element size: 0.005m)

한편, 두 대상물 모두 요소의 크기를 0.005m로 설정하였을 때 최대 변위를 수렴 값으로 간주할 수 있었는데, Table 2에서 확인할 수 있듯 CR 기법으로 해석한 경우가 NASTRAN 대비 요소 수가 적을 때 수렴 값에 빠르게 수렴하는 것을 확인할 수 있었다.

Table 2.

Comparison of the maximum displacement

Mesh Element Size(m) CR-Shell NASTRAN
0.05 2.9043×10-3 (97.17%) 2.5161×10-3 (83.28%)
0.025 2.9583×10-3 (98.98%) 2.8182×10-3 (93.27%)
0.01 2.9825×10-3 (99.79%) 2.9916×10-3 (99.01%)
0.005 2.9888×10-3 3.0124×10-3

2.2 해석 대상

이 장에서는 본 연구에서 사용될 모델의 생성과 해석에 필요한 물성값과 해석 조건을 설명하였다.

2.2.1 모델링과 경계 조건

본 연구에서는 실제 파손이 일어났던 야마하 사의 FJR-1300 프레임의 형상을 참고하여 모델링하였다. 모터사이클 프레임이 경량화를 위해 속이 빈 쉘 프레임을 이용하였으므로 쉘 요소로 모델링하였으며, 경계 조건은 Fig. 10의 좌측 하단 4지점이 고정단으로 변위와 회전이 구속되었다. 한편 수직 방향 하중이 후방프레임의 탑케이스 적재 부분인 화살표가 있는 곳에 인가되었다. 하중은 제조사에서 판매하는 탑케이스의 최대 용적 중량+자체 중량을 합하여 합계 5.5kgf로 설정하였다.

http://static.apub.kr/journalsite/sites/jcoseik/2020-033-04/N0040330401/images/Figure_jcoseik_33_04_01_F10.jpg
Fig. 10.

Main geometry

2.2.2 물성값

재질은 실제 모터사이클 프레임에서 가장 많이 사용되는 Al6061-T6를 사용한 것으로 가정하였으며, 등방성으로 가정하였다. 이 때 물성값은 Table 1과 같다.

2.2.3 설계변수해석 모델 및 조건

설계변수는 응력 집중 부위의 보강을 위하여 통상적으로 많이 사용되는 두께를 늘리는 방안과 보강용 자재를 도입하는 방안을 고려하여 Table 3과 같은 설계변수해석 조건을 설정하였다. 이때, 증가한 중량은 일정하게 유지하기 위해 증가한 체적이 일정하게 유지되도록 각 형상을 설정하였으며, 보강용 자재는 Fig. 11과 같이 설정되었다. 하중은 비교의 편의성을 위해 500N으로 설정하였다.

Table 3.

Case models for parametric analysis

Case 1 Thickness change of weak plate (2.0mm → 3.0mm)
Case 2 A gusset install (Thickness 8mm)
Case 3 Case 1+Case 2 (Plate Thickness 2.0mm → 2.5mm, gusset thickness 4mm)
http://static.apub.kr/journalsite/sites/jcoseik/2020-033-04/N0040330401/images/Figure_jcoseik_33_04_01_F11.jpg
Fig. 11.

Gusset

2.3 연구 진행

이 장에서는 본 논문에서 수행한 연구 흐름도를 제시하고, 진행 과정 및 결과를 설명한다. 연구는 프레임 형상의 정적/동적 해석 부분, 설계변수해석 부분으로 구성되어 있다.

Fig. 12에 나타난 것과 같이 각 해석 부분들은 해석 입력 파일 작성과 정적 해석, 그리고 동적 해석의 단계를 거치며 그 상세한 방법들은 다음과 같다.

http://static.apub.kr/journalsite/sites/jcoseik/2020-033-04/N0040330401/images/Figure_jcoseik_33_04_01_F12.jpg
Fig. 12.

Research flow chart(main)

2.3.1 해석 입력 파일 작성

모든 형상의 전처리 과정은 SolidWorks로 만든 모델을 EDISON에서 제공하는 CSD Pre 전처리기를 활용하여서 수행하였다. Fig. 13과 같이 쉘 요소로 이산화를 수행하였으며 경계 조건은 2.2절에서 설명한 바와 같다.

http://static.apub.kr/journalsite/sites/jcoseik/2020-033-04/N0040330401/images/Figure_jcoseik_33_04_01_F13.jpg
Fig. 13.

Discretized geometry of frame model obtained by CSD pre

http://static.apub.kr/journalsite/sites/jcoseik/2020-033-04/N0040330401/images/Figure_jcoseik_33_04_01_F14.jpg
Fig. 14.

Stress distribution by CR shell static analysis

2.3.2 정적 해석

정적 해석은 EDISON을 통한 CR 쉘 정적 해석기법으로 수행되었다. 변위 분포를 확인한 후, ANSYS 정적 구조해석 모듈과 비교하여 검증을 완료하였다. 이후 변위 분포를 통해 응력을 계산하는 CR 쉘 응력 계산 기법을 이용하여 응력 분포를 확인하였다.

2.3.3 동적 해석

동적해석의 경우, EDISON을 통한 CR 동적 해석기법을 이용하였다. 구조물에 수직 방향의 step function을 가하여 변위 응답을 계산하였다. 이후 MATLAB으로 변위 응답을 FFT를 통해 분석하였고, 고유진동수를 예측했다.

2.4 연구 결과

2.4.1 정적 해석

정적 해석결과 응력의 강도가 높게 측정된 부위들이 실제 파손이 발생하는 부위와 일치했다. 그러나 최대 응력은 2.4×107Pa로 파손이 발생하기에 충분한 응력이 발생하지 않았다.

2.4.2 동적 해석

수직 방향으로 step function의 하중을 가한 후, Co-rotational 기법으로 쉘 요소 동적 해석을 통해 시간에 따른 응답을 분석하였다. FFT를 거쳐 주파서 대역으로 응답을 분석한 결과는 Fig. 15와 같다. 그래프에서 보듯 수직 방향 거동에 대해 1차 고유진동수가 78Hz로 계산되었다. 이때 모드 형상은 Figs. 16~17과 같이 표현되었다.

http://static.apub.kr/journalsite/sites/jcoseik/2020-033-04/N0040330401/images/Figure_jcoseik_33_04_01_F15.jpg
Fig. 15.

Vertical motion response

http://static.apub.kr/journalsite/sites/jcoseik/2020-033-04/N0040330401/images/Figure_jcoseik_33_04_01_F16.jpg
Fig. 16.

1st mode shape

http://static.apub.kr/journalsite/sites/jcoseik/2020-033-04/N0040330401/images/Figure_jcoseik_33_04_01_F17.jpg
Fig. 17.

1st mode shape

2.4.3 설계변수 해석

2.2.3에서 제시했던 세 가지 형상에 대해 정적 해석과 동적 해석을 수행하였다. Figs. 18~21에서 알 수 있듯, 최대 응력의 크기는 Case 3> Case 2> Case 1의 순서로 나타났다.

http://static.apub.kr/journalsite/sites/jcoseik/2020-033-04/N0040330401/images/Figure_jcoseik_33_04_01_F18.jpg
Fig. 18.

The stress distribution of the original model under 500N load

http://static.apub.kr/journalsite/sites/jcoseik/2020-033-04/N0040330401/images/Figure_jcoseik_33_04_01_F19.jpg
Fig. 19.

The stress distribution of the 1st case model

http://static.apub.kr/journalsite/sites/jcoseik/2020-033-04/N0040330401/images/Figure_jcoseik_33_04_01_F20.jpg
Fig. 20.

The stress distribution of the 2nd case model

http://static.apub.kr/journalsite/sites/jcoseik/2020-033-04/N0040330401/images/Figure_jcoseik_33_04_01_F21.jpg
Fig. 21.

The stress distribution of the 3rd case model

한편, Figs. 22~24에서 보이듯 각 Case에 따라 응력이 집중되는 부분들이 달라졌다. 판의 두께를 증가시킨 경우에는 실제 파손 케이스와 동일하게 판의 모서리 부분에 최대 응력을 확인할 수 있었으며, 보강용 자재를 장착한 경우 모서리의 응력은 감소하였으나 보강용 자재가 장착된 부분에 높은 응력이 확인되었다. 이는 보강용 자재가 모서리 부분의 강성을 높여주나 보강용자재가 고정된 다른 부분에 그 하중을 전달하기 때문으로 보인다. 한편, 고유진동수의 경우 원 형상과 동일한 결과가 얻어졌다. 이는 전체 프레임에 비해 상대적으로 작은 설계 변경이 있었기 때문으로 보인다.

http://static.apub.kr/journalsite/sites/jcoseik/2020-033-04/N0040330401/images/Figure_jcoseik_33_04_01_F22.jpg
Fig. 22.

Maximum stress point in the 1st case model

http://static.apub.kr/journalsite/sites/jcoseik/2020-033-04/N0040330401/images/Figure_jcoseik_33_04_01_F23.jpg
Fig. 23.

Maximum stress point in the 2nd case model

http://static.apub.kr/journalsite/sites/jcoseik/2020-033-04/N0040330401/images/Figure_jcoseik_33_04_01_F24.jpg
Fig. 24.

Maximum stress point in the 3rd case model

3. 결 론

본 연구에서는 CR shell analysis를 이용하여 단순화된 프레임 형상의 3차원 쉘 구조에 대해 요소개수에 따른 계산 결과를 MSC NASTRAN과 비교하였다. 또한, 3차원 쉘 구조로 모델링 된 모터사이클의 후방프레임에 대한 정적/동적 구조해석을 진행하고 상용 소프트웨어인 ANSYS와의 검증을 완료하였다. 이후 2가지 설계변수가 적용된 3가지 형상에 대해 해석을 진행하고 결과를 비교하였다. 이에 대한 연구 결과는 다음과 같다.

1)단순화된 프레임 형상을 통해 비교해 본 결과 CR 기법은 NASTRAN과 비슷한 결과를 보여 주면서도 상대적으로 적은 요소 수에서 높은 정확도를 보여 주며 프레임 해석에서 더욱 경제적인 방법임을 입증하였다.

2)정적 해석결과 해당 프레임 형상에서 실제 사례에서 관찰되었던 파손 부위에 최대 응력 2.4×107Pa이 확인되며 집중됨을 확인할 수 있었다.

3)모드 해석을 통해 하중 조건에서 78Hz의 부근에서 공진을 일으킬 수 있음을 확인하였다. 이는 제조사에서 권고한 130km/h 부근에서의 엔진 회전속도(5단, 4,400RPM)과 근접한다. 따라서, 주행 중 프레임의 부분적인 공진 발생으로 인한 수명 저하 가능성을 고려할 수 있다.

4)설계변수 해석결과, 두께를 증가시키는 것과 보강용 자재를 설치하는 것은 최대 응력을 낮춰주며, 특히 보강용 자재는 최대 응력이 발생하는 부위를 이동시키는 결과를 확인할 수 있었다. 한편 최대 응력감소 효과는 두께를 증가시키는 것과 보강용 자재를 설치하는 것을 동시에 적용했던 Case 3이었으며, 이를 통해 최대 응력이 Case 1과 Case 2에서 발견되었던 두 부분으로 분산되며 낮아짐을 확인할 수 있었다.

따라서 본 연구의 결과를 종합하면 프레임 해석에 있어 CR 쉘 요소를 활용한 기법이 효과적이며, 모터사이클 후방프레임에서 특정 부위에 응력이 집중되고 실제 주행 시 엔진에서 가해지는 진동수가 프레임의 고유진동수와 근접함을 알 수 있었다. 그러나 하중만으로는 파괴에 이르기 충분치 않으므로 실제 파괴 원인을 규명하기 위해서는 공진이 미치는 영향을 고려해야 함을 추론할 수 있었다. 또한, 설계변수해석을 통해 두께를 증가시키는 것보다는 보강용 자재를 설치하는 것이 동일한 중량 증가 시 응력감소에 효과적임을 확인할 수 있었으며, 보강용 자재만 설치하기보다는 두께 증가와 함께 보강용 자재를 설치하는 것이 더욱 효과적임을 확인할 수 있었다. 본 연구에서는 요소를 절감할 수 있는 CR기법과 더불어, 클라우드 서비스로 이루어지는 EDISON 소프트웨어를 사용하여 동시에 여러 계산을 수행할 수 있어 시간을 절감할 수 있었다. 본 연구의 결과를 종합하여 모터사이클 후방프레임의 설계에 응용한다면 고속 주행에서 프레임 파손을 피할 수 있는 설계를 고안하는데 유용할 것으로 생각되며, 후속 연구에서는 본 연구에서 확인된 판의 두께 증가와 보강용 자재 설치가 응력감소에 미치는 영향을 활용하여 프레임의 응력감소를 극대화하는 최적화 연구를 수행하고자 한다.

Acknowledgements

본 연구는 정부(과학기술정보통신부)의 재원으로 한국연구재단 첨단 사이언스·교육 허브 개발 사업(NRF-2011-0020576)의 지원을 받아 수행되었습니다.

References

1

Cho, H., Shin, S. (2015a) Three-dimensional Precision Analysis for Flapping Wing Structure using CR Shell Elements, Annu. Conf. Korean Soc. Aeronaut. & Space Sci., pp.942~945.

2

Cho, H., Shin, S. (2015b) Development of Nonlinear Triangular Plane Element Based on Co-rotational Framework, J. Comput. Struct. Eng. Inst. Korea, 28(5), pp.485~490.

10.7734/COSEIK.2015.28.5.485
3

Joo, H., Cho, H., Kim, S., Shin, S. (2017) Large-Sized Three- Dimensional Structural Analysis using OPT-DKT Shell Element based on FETI-Local Method, Annu. Conf. Comput. Struct. Eng. Inst. Korea.

4

Kim, T. (2018) Yamaha Motorcycle Frame Damaged..."Defect" VS "Fault", MBN News.

5

Lee, Y.W., Ha, S.Y., Kwon, J.H. (2016) Parametric Analysis for Structural Stiffness Enhancement of Motorcycle Framel, J. Trans. Korean Soc. Automot. Eng., 24(5), pp.612~617.

10.7467/KSAE.2016.24.5.612
6

Park, B. (1998) Vibration Characteristics of a Motorcycle Body, J. Korean Soc.r Precis. Eng., 15(1), pp.169~176.

페이지 상단으로 이동하기