전산소성역학 기말 프로젝트

2022-21559 곽민준

main program

1. 각 loading mode에 대해서 strain 정보 읽음

   1. LE11, LE22, LE12

2. elastic property 기입

   1. E, v

3. hardening property 기입

   1. f=sigma_e-sigma_0 에서, sigma_0 이 equivalent plastic strain(EPS)의 함수
   2. 이 프로그램에서는 linear hardening가정
   3. sig_Y, k: initial yield stress, coeff. of linear hardening

4. state variable initiating

   1. EPS: equivalent plastic strain

5. stress initiating

   1. plane stress인 stress=0벡터로 만듦

6. loop over the strain

   1. 이전 step strain과 현 step strain으로 행렬 생성하여, 전산소성역학 기말 프로젝트에 대입
      → stress updated, EPS, sig_0 출력

   2. S11, S12, S22에 저장하여 다시 다음 strain으로 넘어감

7. 모든 step 끝나면 그래프 출력

전산소성역학 기말 프로젝트

1. 입력값
   1. E: 200Gpa, v: 0.3
   2. sig_Y: 200Mpa, k: 2000Mpa
   3. stress
   4. strain_prev, strain_curr
   5. eps
   6. tol1, tol2 = 1E-8, 1E-4
2. E, v 값을 전산소성역학 기말 프로젝트에 대입하여 elastic stiffness matrix (C_ela) 반환
3. 탄성예측 시도

   1. 저장된 stress에 현 step과 이전 step 사이 delta strain을 이용해 완전 탄성을 가정한 탄성 예측 stress 계산
      stress = stress + C_ela*dstrain;

      $\Delta {\mathbf{\sigma }}_e=\mathbf{C}\Delta {\mathbf{\epsilon }}_t$

   2. stress0 (탄성예측된 stress값)으로 저장해둠
      stress0 = stress

4. 전산소성역학 기말 프로젝트 함수에 탄성예측된 stress를 대입하여 equivalent stress(sig_e) 반환
5. equivalent strain을 전산소성역학 기말 프로젝트 함수에 hardening parameter들과 같이 대입하여 yield stress를 반환
6. 이 equivalent stress(sig_e)값이 yield stress(sig_0)보다 작을 경우 elastic 구간이므로 별도의 계산 없이 다음 루프로 넘어가지만 더 클 경우 plastic region에 들어선 것으로, 여기부턴 뉴턴랩슨 기반의 소성 보정이 필요

   1. f = sig_e - sig_0, f>0이면 소성 보정

      $f=({\sigma }_x^2+{\sigma }_y^2-{\sigma }_x{\sigma }_y+3{\tau }_{xy}^2{)}^{1/2}-{\sigma }_0$

   2. f<0이면 N-R 없이 바로 다음 strain으로 loop가 넘어간다.

7. 소성 보정을 위한 전산소성역학 기말 프로젝트사용(Euler backward scheme)

   1. iteration을 통해 delta lambda(plastic multiplier) 값을 보정

   2. d(delta lambda) 값을 더해가면서 dlam값 update, stress도 마찬가지로 update

   3. 만약 residual값과 f값이 tolerance보다 작으면 yield function 위에 있다고 판단, loop 종료

      $$\begin{gathered} \text{r}={\mathbf{\sigma }}_C-({\mathbf{\sigma }}_B-\Delta \lambda \text{C}{\text{a}}_c) \\ |\mathbf{r}{|}^2<tolerance1 \end{gathered}$$

8. 소성 보정이 완료된 delta lambda값을 eps에 update
9. 다음 step으로 넘어감

전산소성역학 기말 프로젝트

1. 스칼라값 E, v값을 받아 elastic stiffness matrix(C_ela) 반환

$C=\frac{E}{1-{\nu }^2}\left[\begin{array}{c}1\nu 0\\ \nu 10\\ 00\frac{1-\nu }{2}\end{array}\right]$

전산소성역학 기말 프로젝트

1. stress [sigma_11;sigma_22;sigma_12]를 받아서 equivalent stress 반환

   ${\sigma }_v=\sqrt{{\sigma }_x^2-{\sigma }_x{\sigma }_y+{\sigma }_y^2+3{\tau }_{xy}^2}$

2. a는 yield function 위에서의 normal 방향이며, gradient of yield function

   $$\begin{gathered} f=\sqrt{{\sigma }_x^2-{\sigma }_x{\sigma }_y+{\sigma }_y^2+3{\tau }_{xy}^2}-{\sigma }_0={\sigma }_e-{\sigma }_0 \\ \Delta {\mathbf{\epsilon }}_p=\Delta \lambda (\frac{\partial f}{\partial \mathbf{\sigma }})=\Delta \lambda \mathbf{a} \\ \Rightarrow \left(\begin{array}{c}{\dot{\epsilon }}_{px}\\ {\dot{\epsilon }}_{py}\\ {\dot{\epsilon }}_{pxy}\end{array}\right)=\frac{\dot{\lambda }}{2{\sigma }_e}\left(\begin{array}{c}2{\sigma }_x-{\sigma }_y\\ 2{\sigma }_y-{\sigma }_x\\ 6{\tau }_{xy}\end{array}\right) \end{gathered}$$

전산소성역학 기말 프로젝트

1. equivalent plastic strain(eps), coeff. of linear hardening(k), initial yield stress(sig_Y)를 받아 yield stress(sig_0) 반환
2. 여기선 linear hardening을 사용

${\sigma }_0({\epsilon }_{ps})={\sigma }_Y+k{\epsilon }_{ps}$

전산소성역학 기말 프로젝트

1. TODO: delta(delta lambda), delta stress 계산

2. 4 Unknowns: delta lambda(1 scalar unknown), stress at C(3 unknown values on the yield surface)

   $\{\Delta \lambda ,{\mathbf{\sigma }}_C\}where{\mathbf{a}}_C=[\partial f/\partial \mathbf{\sigma }{]}_C$

3. 4 equations

   1. r vector에서 3equation, f값에서 1 equation
      강의자료의 경우 3D(7x7)였으나 2D(4x4)로 수정하여 사용하였음

      $0=\left(\begin{array}{c}{r}_1\\ r_2\\ r_3\\ f\end{array}\right)$

      $\{\begin{array}{l}\mathbf{r}={\mathbf{\sigma }}_C-({\mathbf{\sigma }}_B-\Delta \lambda \mathbf{C}{\mathbf{a}}_C)=0\\ f({\mathbf{\sigma }}_C)={\sigma }_e({\mathbf{\sigma }}_C)-{\sigma }_0({\epsilon }_{ps})=0\end{array}$

   2. simplified matrix equation

      $( \mathbf{r}=\boldsymbol{\sigma}_{C}-\left(\boldsymbol{\sigma}_{B}-\Delta \lambda \mathbf{C} \mathbf{a}_{C}\right)=0)$

4. matrix 각 성분 계산

   1. upper left

      1. r에서 trial stress(sigma_B)는 fixed value이므로,

         

   2. upper right

      $\frac{\partial \mathbf{r}}{\partial \Delta \lambda }=\mathbf{Ca}$

   3. lower left

      1. a의 definition에서,

      ${\mathbf{a}}_C=[\partial f/\partial \mathbf{\sigma }{]}_C\Rightarrow \frac{\partial f}{\partial \sigma }=\mathbf{a}$

   4. lower right

      1. Von Mises에서 delta lambda = epsilon_ps 이므로,

         $\frac{\partial f}{\partial \Delta \lambda }=\frac{\partial {\sigma }_0}{\partial {\epsilon }_{ps}}=A^{\prime }$

      2. A’값은 이 코드에선 k값으로 사용됨

5. update를 위한 delta 값들 계산

   1. 단순화를 위해, Q matrix를 정의

      $\mathbf{Q}\doteq \mathbf{I}+\Delta \lambda \mathbf{C}\left(\frac{{\partial }^{2}f}{\partial {\mathbf{\sigma }}^2}\right)$

   2. 아래 equation을 풀어 delta값 계산

      0=\left(\begin{array}{l}\mathbf{r} \\ f\end{array}\right) +\left(\begin{array}{l}\Delta \mathbf{r} \\ \Delta f\end{array}\right)\\ \\ 0=\left(\begin{array}{l}\mathbf{r} \\ f\end{array}\right)+\left[\begin{array}{cc}\frac{\partial \mathbf{r}}{\partial \boldsymbol{\sigma}} & \frac{\partial \mathbf{r}}{\partial \Delta \lambda} \\ \left(\frac{\partial f}{\partial \boldsymbol{\sigma}}\right)^{T} & \frac{\partial f}{\partial \Delta \lambda}\end{array}\right]\left(\begin{array}{c}\Delta \boldsymbol{\sigma} \\ \Delta^{2} \lambda\end{array}\right)\\ \left[\begin{array}{cc}

6. 계산된 delta 값을 update 후 f, r값 tolerance와 비교하여 수렴 여부를 파악

결과 분석

※ 각 case에서 왼쪽 위에서부터 각각 11, 22, 12 방향에 대하여 stress를 코드를 통한 계산한 값을 y, 아바쿠스로 추출된 데이터를 x로 하여 plot한 것임

1. plots

   1. Uniaxial

      

   2. Biaxial

      

   3. Shear

      

2. 세 케이스 모두 정상적으로 수렴하였으며, 데이터를 보았을 때 모든 mode에서 코드를 통해 도출된 stress값과 레퍼런스 아바쿠스 stress 값이 대체로 일치하는 것을 확인할 수 있다. uniaxial의 경우 S11에서 약간 값이 튀긴 했으나 거의 0에 가까운 값이 도출되었으며, 이는 아바쿠스 데이터와도 거의 유사하다.
3. iteration의 경우 total iteration update 코드를 추가하여 확인한 결과 uniaxial에서 총 555번, biaxial의 경우 600번, shear의 경우 550번으로 집계되었다.
4. 효율 평가를 위해 strain increment를 현재 값에서 2배로 하여 계산했을 때 total iteration 횟수는 더 줄어든 318번으로 집계되었다.