2018년 1월 2일 화요일

Stochastic Mechanics


지반공학을 위한 수치해석, 정말 쉽게 한번 접근해 보자.


시작하긴 겁나지만 그 막강함에 안할수가 없는 Stochastics에 대해 요점정리 해본다.



Fundamental concepts: 기본 개념


Stochastic? Monte Carlo? Markov-chain? Bayesian? Probabilistic? Uncertainty? Reliability?????

천천히 시작해보기 위해 wikipedia에서 정의부터 확인해 보자.


Definition

Stochastics =  Theory of statistical inference

: 추계학 = 추측 통계학 =  "모집단에서 임의로 추출한 표본에 따라 모집단의 상태를 추측하는 학문"

어렵다. 정작 필요한 Uncertainty 혹은 Reliability의 개념을 정의로부터 찾을 수 없기 때문에 필자도 Stochastics의 시작이 너무 어려웠다. 하여 다시금 EZ하게 설명해본다.

EZ Definition

Stochastic: Modeling 결과값의 Uncertainty를 구하기 위해 random number를 활용하는 학문

여기서...
Modeling이란 "A model is a theoretical way of understanding a concept or idea."
Uncertainty란 통계학적 정규분포(Normal distribution) 상의 표준편차(Standard deviation)
Reliability란 확률에 입각한 mechanical analysis의 결과이자 최종 목표이다.

정의에서 알 수 있듯이 Stochastic의 핵심은 modeling parameter가 random하다는 것이다. 이 random variable을 백만번 modeling 돌리게 되면 결과값이 정규분포를 띄게 되고 거기서 표준편차값을 구해 Reliability를 계산하는 것이다. 하여, 모델링을 백만번 돌리기 위한 알고리즘이 Markov-chain Monte Calro (MCMC) method이고 이 MCMC에 대한 이해가 정말 중요하다.

Algorithm: MCMC method


MCMC의 이해 또한 wikipedia에서부터 시작해 보자.

"마르코프 연쇄 구성에 기반한 확률 분포로부터 원하는 분포의 정적 분포를 갖는 표본으 추출하는 알고리즘"

이게 뭔 소린가. 차라리 영어로 다시 보자.

"MCMC methods are a class of algorithms for sampling from a probability distribution based on constructing a Markov chain that has the desired distribution as its equilibrium distribution."

즉, 우리가 원하는 Markov chain을 만들기 위해 확률적으로 sampling하는 방법이라는 얘기다. 다시금 EZ하게 설명해보자면

EZ MCMC method

MCMC: 원하는 Markov-chain을 만들기 위한 Monte Carlo, 즉 random number 활용한 반복 연산 방법

여기서...
Monte Carlo method란 random number를 활용한 수천만번의 모집단 생성 방법이고
Markov-chain Monte Carlo란 모집단이 너무 많으니 확률적으로 일부만 추려내자는 것(sampling) 이며
Bayesian probability가 그 추려내는 과정에서의 척도가 되는 것이다.


다시금 정리해 보자면,

Monte Carlo의 randomness와 Bayes' rule의 probability를 결합하여 Markov-chain을 만들어 내는 것이 MCMC method이고, 이 알고리즘을 활용해서 모델링의 Uncertainty를 확인하고 그에 따른 Reliability를 계산하는 것이 Stochastic Mechanics가 된다.


Bayesian Inference

베이지안을 또한 wikipedia에서 검색해보면
"통계적 추론의 한 방법으로, 추론 대상의 사전 확률과 추가적인 정보를 통해 해당 대상의 사후 확률을 추론하는 방법이다." 라고하고 영어로는
"Statistical inference in which Bayes' theorem is used to update the probability for a hypothesis as more evidence or information becomes available."이란다.

한마디로, Bayesian은 조건부 확률을 통한 추측이다.

예를 들자면,
  • 산사태가 일어날 확률이 그냥 추측이고
  • 내일 비가 오는데 산사태가 일어날 확률이 조건부 확률을 통한 추측이다.

그래서 다음의 세가지를 필요로 한다.
  • Prior : 사전 정보, 즉 내일 비가 온다는 정보, 사면 기울기, 토질 저항력 등등
  • Likelihood : 조건식, 즉 토질의 특정 함수비 상태에서의 산사태 확률
  • Posterior : 베이지안의 주 목적 결과값, 즉 비 올 때 산사태 확률

따라서, Posterior는 기존 경험상의 Prior에서 데이터상의 Likelihood를 접목시킨 가장 믿을만한 새로운 기준이 되고, 이 기준에 맞춰서 수천만개의 모집단 중 일부만 추려서 sampling할 수 있게 된다.

Sampling

sampling은 MCMC 알고리즘의 일부 테크닉으로서 여러가지 방법이 있지만, 다음 두가지가 널리 쓰인다.
  • Metropolis-Hastings(MH) sampling algorithm
  • Gibbs sampling algorithm

간단히 얘기해서 MH가 일반적인 알고리즘이고, Gibbs는 특정 조건에 특화된 간단한 알고리즘이다. MH sampling에서는 한가지, probability of acceptance (alpha)를 알아야 한다. 이 alpha는 Bayesian에서 계산되어 나오고, 그 값이 램덤상수 r값을 통해 확률적으로 accept되게 만든 것이다. 물론, 이 accept된 값들은 Markov-chain으로서 저장되어 우리가 원하는 결과값을 만들 수 있게 된다는 말이다.


한 가지 tip: Normal Distribution(정규분포)

알고보니 모든 확률적 방법론들은 모집단이 정규분포임을 가정하고 변수들을 조정한다. 따라서 Stochastics에 대한 개념이해가 됐다면 이미 알고 있다고 간과해버렸던 Normal distribution의 개념 또한 확실히 하고 예제를 다뤄볼 필요가 있다.

정규분포의 가장 중요한 특성은 2개의 매개 변수 (평균, 표준편차)에 의해 모양이 결정되는 연속 확률 분포라는 것이다.  따라서, 수식 몇개만 유념하고 넘어가자.

Notation of normal distribution is:

The probability density of the normal distributions is:

이때, 정규분포가 다차원 공간에 대해 확장된다면,

Notation of multivariate normal distribution is:

The probability density function(PDF) is:

이 PDF가 정말 중요하고 많이 쓰이고 어려워 보이지만 유심히 보면 일반 정규분포 식이랑 거의 똑같다. 당연하다. 다만 분산(variance = standard deviation^2) 이 공분산(covariance) matrix로 바뀌었을 뿐이다.

그럼 covariance matrix는 무엇이냐?

공분산은 2개의 확률변수의 상관정도를 나타내는 값인데 + 혹은 - 값을 가지고 0 이면 두 변수는 독립이라는 말이다. 물론, 상관관계를 나타내기에 부적절해서 -1 과 +1 사이값으로 표준화 시켜놓은 상관계수 (correlation)를 우리는 암암리에 쓰고 있지만 말이다.

쉬운 산수 계산 예를 들자면,
Xi = [0, 1, 2, 3, 4] 이고 Xj = [2, 3, 2, 1, 2] 일 때, 평균은 각각 mi = 2, mj = 2 이다.
이때 Xi의 population 분산 = (4+1+0+1+2)/5 = 2 이고 sample 분산 = (4+1+0+1+2)/4 = 2.5 이다.
공분산 cov(Xi,Xj) = [[-2 -1 0 1 2],[0 1 0 -1 0]] * [[-1 0],[-1 1],[0 0],[1 -1],[2 0]] 로서 사이즈가 (2x5)x(5x2) = (2x2) 매트릭스인 [[2.5 -0.5],[-0.5 0.5]]가 나온다.

그러니까 우선은 이렇게만 알고 넘어가자. ㅡ.ㅡ;;



EZ Example


기본 개념을 이해했다면 MCMC algorithm의 코딩을 위해 Bayes' rule에 입각하여 쉬운 예제를 이해해보자.

To be updated



Probabilistic Inversion


MCMC 실질적 활용은 크게 세가지라고 생각한다.
  • Uncertainty and Reability against structural failure, 왜? 모든 재료는 균질하지 않기 때문에
  • Regression: 회귀분석, 왜냐? 마더네이쳐는 Non-linear 라서
  • Inversion: 역산, 왜냐? 역산은 원래 non-uniquness 라서 ill-posed problem 때문에

이중 필자의 연구와 가장 가까운 Inversion(역산)에 대해 상세하게 다시 설명해보고자 한다. MCMC inversion의 key point는 MCMC 자체의 이해도 이해지만, 본인이 원하는 forward modeling을 MCMC 안에 제대로 접목시키는 것이다.

즉, 코딩을 위한 알고리즘은 실측 데이터(d_obs)의 포맷 확인 후, 그에 맞는 forward modeling code의 테스트로부터 시작되어야 하며, initial 모델링 결과값과의 비교대조 후, MCMC를 적용시켜야 한다.
  1. check observed data (d_obs)
    1. d_obs는 측정을 통한 실측 데이터값이다.
  2. test forward modeling
    1. d = g(m)이 기본 개념으로서 수학적으로 결국 y = f(x)와 같은 개념이다.
    2. 이때 d = d_pred로 모델링을 통한 예상 데이터 값이다.
    3. 결국, 우리는 d_obs와 가장 유사한 d_pred를 만들어서 실제 m값을 추정하는 것이다.
  3. MCMC algorithm
    1. nsampling 값을 정하고 그에 맞춰 모든 변수들의 memory allocation
    2. initial modeling으로부터의 d_pred값을 Markov-chain의 첫 element(m0)로 설정
    3. start constructing candidate model(m')
      1. 이때 보통의 MCMC는 move 혹은 perturb parameter만 이뤄지지만
      2. rj-MCMC의 경우는 birth/death 개념까지 포함된다.
      3. perturbation: m' = m_(k-1)+dm_(k)
    4. compute the probability of acceptance(alpha) of m' from Bayes' rule
      1. likelihood
      2. prior
      3. posterior: P(m_(k)|d)
    5. M-H sampling: random number(r) < alpha 일때만 m'을 accept
      1. r의 적용을 통해서 확률적으로 accept하기 때문에,
      2. local minima에 빠지지 않고 나올 수 있게 된다. (optimization 관련)
      3. 아무튼 m'이 accept되면 m_(k) = m' 하고 k=k+1
    6. inversion의 수렴 확인 후, burn-in point 이후부터 nsampling만큼 Markov-chain 형성
  4. Post processing


친절한 review

통계학과 Bayesian Modeling and Inference 수업의 친절한 syllabus를 통해 우리가 빠트린 부분을 보충해보자.


Bayesian Modeling and Inference - Course topic

  1. Review of distribution theory and basic sampling methods
  2. Fundamentals of Bayesian inference
  3. Conjugate priors in one-parameter families
    1. Conjugate prior 란: posterior distribution 과 same family distribution 을 갖는 prior 이다.
    2. 하지만 이게 왜 중요할까???
  4. Bayesian inference, Credible sets, Marginal likelihood & Bayes factor, Hypothesis testing
  5. Posterior computation with non-conjugate priors
    1. Laplace approximation, Importance sampling, Metropolis-Hastings algorithms
    2. Markov chain Monte Carlo methods, Gibbs & slice sampling, MCMC diagnostics
  6. Bayesian linear & generalized linear model
  7. Bayesian model selection & model averaging
  8. Hierarchical models, Mixure models and dynamic models
    1. Hierarchical modeling 이란: Bayesian method를 사용하여 multiple level로 이루어진 posterior distribution 의 parameters 를 estimate 하는 stastical model 이란다.
    2. 그래서 뭘까???
  9. Decision theory & asymptotics





2015년 11월 19일 목요일

Geotechnical Soil Properties


지반공학을 위한 수치해석, 정말 쉽게 한번 접근해 보자.


지반공학의 절대 기본이 되는 흙의 특성들을 요점정리 해본다.



대전제


: 흙은 Solid (Sand or Clay) + water + air 혼합물로서,


 heterogeneous, anisotropic 한 비선형 거동 물성체이다.



The unified classification method (통일분류법)


Sieve analysis --> S / C 구별 and W / P graded 판단

  • Cu : Coefficient of Uniformity = D60 / D10
  • Cc : Coefficient of Curvature = D30^2 / (D10*D60)
  • if Cu > 6 and 1 < Cc < 3, then SW

Atterberg limit test --> Plasticity Chart
  • LL (wL) : Liquid Limit :차트 가로축
  • PL : Plastic Limit
  • PI (Ip) : Plastic Index = LL - PL : 차트 세로축




etc.
  • e : void ratio = Vv / Vs  <-- 가장 많이 쓰임. 중요.
  • w : water content = Mw / Ms
  • Sr : Degree of Saturation = Vw / Vv
  • n : porosity = Vv / V = e / (1+e)
  • gamma_t : total unit weight = (G+Sr*e) / (1+e) *  gamma_w


Stress & Strain


: 아래에선 기본적인 응력 sigma 를 P 로 표기한다.

Mohr's Circle
  • Shear strength
  • Mohr-Coulomb failure envelope --> c , phi

Effective stress
    • Terzaghi (1936)
      • The effective stress = total stress - pore water pressure
      • P' = P - pwp
    • Earth Pressure at Rest
      • "at Rest" : a stress state of K0 is a constant.
      • K0 ( = P'h / P'v ) : Earth pressure coefficient at rest

    p & q
    • Lambe representation : 2D --> Mohr's Circle 대신 사용
      • p : average stress = 1/2 * (P1 + P3) = p' + pwp
      • q : deviator stress = 1/2 * (P1 - P3) = q'
    • Cambridge representation : 3D --> Triaxial Test 결과플랏 y축에 q 사용
      • p = 1/3 * (P1 + P2 + P3) = p' + pwp
      • q = (P1 - P3) = q'


    Triaxial compression test
    • UU : Unconsolidated Undrained = Quick test
      • 초기 안정만 해석, 보통 시공 직후가 가장 위험하므로
      •     Saturated : c > 0, phi = 0 --> tau : cu = 0.5(s1-s3)
      • Unsaturated : c > 0, phi > 0
    • CU : Consolidated Undrained = Slow test
      • 흙댐, 제방 완공 후 수면 급격히 하강한 경우 안정계산
      • CU: effective stress test --> c' & phi'
    • CD : Consolidated Drained = 어려운 test 라 CU 로 보통 대체
      • 투수계수가 큰 흙에 사용

    Shear strength
    • Clay --> Sensitivity = undisturbed strength / remoulded strenfgh ( if >16 : Quick clay )
    • Dense sand --> Dilatancy
    • Loose sand + saturated + undrained + shear stress --> Liquefaction 

    • 추가 : soil 별 보통의 c 혹은 phi 값
      • Soft clay : c = 20 ~ 40 kPa
      • Firm clay : c = 50 ~ 75 kPa
      • Loose sand : phi = 27 `
      • Dense sand : phi = 35 `



    Seepage (물의 침투)


    Darcy's Law : Q = k i A


    • k : coefficient of permeability = hydraulic conductivity
      • 1st k 유추식 (Hazen's eq) : k = 0.01 * (D10)^2 [m/s]
      • 2nd k 유추식 (Kozeny-Carman's eq) : k = f(phi) * d^2 * f(s)
      • Flow net : Q = k * h  * Nf / Nd  --> estimate quick conditions (boiling & piping)
    • Assumptions
      • laminar flow
      • isotropic, homogeneous
      • rigid soil structure
      • saturated


    Hydraulic test

    • Constant head test : Q / A = k * L/dH
    • Falling head test : Q dt = - a dh = k * h/L * A dt

    etc.
    • Bernoilli's eq : h = z + u/gamma_w + v^2/(2g)
    • Laplace eq : d2h/dx2 + d2h/dy2 = 0




    Compaction (다짐)


    : 공기 빼기. 단, 100% 제거는 불가능.

     최적의 다짐을 위해선 적정한 물 (water film)이 필요함.
    • Dry density (건조중량) 식이 중요
      • rho_d
      • = rho / (1 + w)
      • = Gs * rho_w / (1 + e)
      • = Gs * (1 - Ar) * rho_w / (1 + w * Gs)
    • Find e, n ( = e/(1+e) )
    • Degree of Saturation 식도 중요
      • Sr 
      • = w * Gs / e
    • RC : Relative Compaction = rho_d (field) / rho_d_max (lab test)



    Consolidation (압밀) : S = Se + Sc + Ss


    : 물 빼기. 문제는 Clay 의 low permeability

    • Clay --> OCR : OverConsolidation Ratio = Past Sv / Present Sv > 1
    •         --> sensitivity


    Initial (elastic) settlement (탄성 침하)


    • Elastic half-space method
      • Assumes homogeneous, elastic, and isotropic material goes on to infinity below a plane surface.
    • Boussinesq
      • Point loads : Pz = Q / z^2 * Ip
      • Line loads --> bulb of pressure until z = about 3B
      • Strip loads
      • Area loads
    • The Newmark influence chart
      • The vertical stress point is required at the center of the chart.
      • Pz = 0.005 * N * q
    • Schleicher (1926)
      • Elastic (initial) settlement : Se = dP * B * (1-v^2) / Es * Ip
      • Se is direct proportion to both dP and B.
    • 추가 : soil 별 보통의 Es 와 v 값
      • Soft clay : Es = 1.8 ~ 3.5 MPa, v = 0.15 ~ 0.25.
      • Hard clay : Es = 6 ~ 14 MPa, v = 0.2 ~ 0.5.
      • Loose sand : Es = 10 ~ 28 MPa, v = 0.2 ~ 0.4.
      • Dense sand : Es = 35 ~ 70 MPa, v = 0.3 ~ 0.45.
    • 추가 : Schmertman & Hartman method (1978)
      • Elastic (initial) settlement : Se = C1 * C2 * q * dz * (Iz / Es)
      • C1 : correction factor for depth of foundation embedment
      • C2 : correction factor to account for creep in soil


    Primary consolidation settlement (압밀 침하)


    • 우선 로그스케일 플랏 [ (e) vs. log(P) ] 를 그린다.
    • OCR = Pc / Po > 1
      • Pc : 선행압밀하중 = 과거 최대 하중
    • Cc : Compression index : 압축지수 
      • Cc = 기울기 of log scale
      • Cc = -(e1 - e2) / (logP1 - logP2)
      • 참고) Skempton 경험식 : Cc = 0.007 * (WL - 10) 
    • Cr : Re-compression index
    • Settlement (Sc) = Cc / (1+e) * H * log[ (P0 + dP) / P0 ]
      • if P0 < Pc < P0 + dP
      • then Sc = Cs/(1+e)*H*log(Pc/P0) + Cc/(1+e)*H*log((P0+dP)/Pc)

    • 다른 방법도 있다.
    • av : coefficient of compressibility : 압축계수 = 기울기 of natural scale
      • av = -(e1-e2) / (P1-P2)
      • mv : coefficient of vomule change: 체적변화율 = av / (1+e0)
    • Settelement (Sc) = mv * del(P) * H

    • 시간에 따른 압밀 침하량 변화 플랏 [ dH vs Time ] 을 그린다.
    • Cv : Coefficient of (Vertical) consolidation : 압밀계수
      • from Terzaghi 1D eqution : du / dt = Cv * d2u / dz2
      • 두 가지 방법
        • root (t) : Taylor : Cv = T90 * H^2 / t90 = 0.848 * H^2 / t90
        • log (t) : Casagrande : Cv = T50 * H^2 / t50 = 0.197 * H^2 / t50
    • U : Degree of consolidation
      • U vs Tv 그래프를 이용
      • 간단히 계산
        • for U < 0.6 --> Tv = pi/4 * U^2
        • for U > 0.6 --> Tv = -0.933 * log(1 - U) - 0.085
    • Consolidtion time (t) = Tv * H^2 / Cv


    • 예를 들어보자. <--- 기억이 안남
      • mv,.... 등이 주어진다.
      • 총 침하량 Sc 를 구한다.
      • 허용 침하 S 와 구해진 Sc 로 U 를 구한다.
      • U 로 Tv 를 구하고 Cv 를 구해서,
      • 최종 시간을 구했었는데...

    Secondary compression settlement (2차 압밀) = Creep



    • Because of the plastic adjustment of soil fabrics
    • Ca : Secondary compression index
      • Ca = - de / log(t1 - t2)
    • Settlement (Ss) = Ca * H * log(t2 / t1)













    Geomechanics



    지반공학을 위한 수치해석, 정말 쉽게 한번 접근해 보자.

    한번은 들어봤지만 한번에 알수가 없는 지반탄성체에 대해 요점정리 해본다.



    Mathematical Preliminaries: 기본 수학


    Tensor (in Euclidian Space)

    • ( a,b ) = scalar                   = integer, potential
    • ( a,b ) = 1st order tensor    = 3x1 matrix vector :  3 components
    • ( A,B ) = 2nd order tensor  = 3x3 matrix            :  9 components
    • ( A,B ) = 3rd order tensor  = 3x3x3 matrix        : 27 components
    • A,B ) = 4rd order tensor  = 3x3x3x3 matrix    : 81 components

    Product
    • inner product                     = · v = c     : projection, eg) (1,0,0) · (0,1,0) = 0+0+0
    • vector (cross) product       = u x v = w   : area,          eg) (1,0,0) x (0,1,0) = (0,0,1)
    • triple scalar (box) product = (u x v) · w  : volume
    • tensor product (dyad)       = u v = U  or  A ⊗ B = D
    • trace                                 = tr(U) = tr(u ⊗ v) = · v
    • dot product                       = A · B = C = Aik Bkj
    • double contraction            = A : B = D = Aij Bij

    Del: Nabla operator (∇)
    • gradient:      ∇ p     = u = ( dp/dx, dp/dy, dp/dz )
    • divergence: ∇ · u   = a = du/dx + du/dy + du/dz
    • curl:             ∇ x u  = v = [ (du3/dy - du2/dz), (du1/dz - du3/dx), (du2/dx - du1/dy) ]
    • Laplacian:    ·∇ p  = d2p/dx2 + d2p/dy2 + d2p/dz2
      • Divergence of Gradient (p), 즉 2차 미분된 scalar 값
    • Hessian:     H(p) = J(u) = J(∇ p)
      • Jacobian of Gradient (p), 이건 2차 미분된 tensor 값


    Divergence theorem (Gauss's theorem): 부피 변화는 전체 면적 flux 변화량 합이다.





    Stoke's theorem: Cauchy's fomula 의 근간. 점에서의 값은 그 점 주변 둘러싼 닫힌 곡선 위에서 적분한 것과 같다??




    Kinematic


    가장 기본은 운동 (혹은 변형) 의 기술방법을 구별해 생각하는 일.
    • Lagrangian Description (X) : Material coordinates: 물체 따라 같이 움직이는 기술방법
    • Eulerian Description (x) : Spatial coordinates: 고정된 관찰자로서의 기술방법

    Deformation gradient tensor (F)
    • dx = F dX
    • F = R U = Rotation tensor x Stretch tensor (symmetric)
    • C = F' F : Right Cauchy-Green = Material deformation tensor
    • B = F F' : Left Cauchy-Green = Finger tensor = Spatial deformation tensor
    • E = 1/2 (C-I)       : Green (Lagrangian) strain tensor
    • e = 1/2 (I - invB) : Almansi (Eulerian) strain tensor
    • 하지만 결국 Cauchy's infinitesimal strain tensor = 1/2 * (u_j,i + u_i,j)

    이게 왜 필요한가?


    구조물 모니터링에서 우리가 얻을 수 있는것은 시간에 따른 좌표값 뿐.
    그에 따른 displacement 와 strain 을 갖고 stress state 를 제대로 판단해야 파괴 예측 가능.




    Stress


    Cauchy stresses
    • T(n) = s_ji * n_j = (3x1) = (3x3)*(3x1)
      • T(n) : traction vector
      • s : Cauchy (true) stress tensor
      • angle between T(n) & n_j = inv( cos( T· n / |T||n| ) )
    • Moment of equilibrium (about x3-axis) : s12 = s21 --> symmetric


    Principal stresses (주응력)
    • Mohr's circle
    • Eigenvalue : det (s -lamda * d_ij) = 0 <--  어렵게만 배우던 고유값의 실제 사용예
    • Invariants (I1, I2, I3) : - s^3 + I1*s^2 - I2*s + I3 = 0

    Stress-Deviation tensor
    • stress tensor = hydrostatic stress tensor (s_m) + deviator stress tensor (s')
    • s'_ij = s_ij - s_m * I
      • s_m = 1/3 * trace(s)
      • I = identity matrix
    • It is very important in describing the plastic behavior.
    • 주의 : triaxial test 에서 보통 쓰던 deviator stress = s1 - s3 랑은 또 조금 다름.
    • 주의 : s' 의 Invariants = (J1, J2, J3) 로 (I1, I2, I3) 와 다름. J1 = 항상 0.

    Constitutive Equations (구성방정식) : s_ij = C_ijkl * e_kl
    • C_ijkl = 4th order tensor : 3^4 = 81 elements
      • s = C e, where C : stiffness tensor
      • e = S s, where S : compliance tensor
    • Symmetrization : 81 - 3*9 - 3*6 = 36 elements
      • Thus, anisotrophy constitutive equation : [s] (6x1) = [C] (6x6) * [e] (6x1)
    • Isotrophiy constitutive equation : s_ij = lamda * e_kk * d_ij + 2 * mu * e_ij
      • s_11 = lamda*(e_11+e_22+e_22) + 2mu*e_11
      • Thus, [C] (6x6) = [ [lamda+2mu, lamda, lamda, 0, 0, 0] ...
    • These can be transformed with ( K, E, v, G (=mu) ) in engineering.
    • 추가 : Compatibility conditions = No gaps or overlaps during deformation.
    • 즉, 적합방정식은 복잡해보여도 결국은 3차원 Strain 들의 관계식.

    한마디로 결론은?


    수치해석에서 등방성 탄성체를 가정했다면, Hooke's Law 는 딱 두 문자로 표현된다.
    바로 Lame's parameters : lamda 와 mu .
    복잡하게 Young's modulus (E), Poisson's ratio (v), bulk modulus (K), shear modulus (G) 로도 표현되곤 하는데, 결국은 똑같은거다. 더 배우고 싶다면 구조방에서...... ㅡ.ㅡㅋ




    Failure Criteria (4 classic models)



    Brittle material : Elasticity
    • 1st. Mohr-Coulomb theory : t = c + s * tan(phi) 
      • t : shear strength
      • c : cohesion
      • s : normal stress
      • phi : angle of internal friction.

    Ductile material : Plasticity

    Beyond the elastic range yield occurs, Elastic and Plastic deformations are assumed to happen concurrently (Total strain e = eElastic + ePlastic) . Plasticity describes the deformation of a material undergoing non-reversible changes of shape in response to applied forces.



    소성변형 들어가기에 앞서서..



    1. EPP model : Elastic Perfectly Plastic model = no work hardening.
    2. 그간 써오던 total accumulated strain 대신, increments in strain 을 사용한다. (loading, unloading 으로 인한 stress-strain 관계가 더이상 1:1 대응이 안되기 때문)
    3. Principal plastic strain increments (d ePlastic) 와 principal deviator stress (s') 는 비례
    4. Flow rule : plastic stress-strain law (between d ePlastic & Cauchy stress (s))
    5. Prandtl-Reuss equations : full elastic-plastic stress-strain relations (w/ Hooke's law)



    • 2nd. Tresca yield criterion
      • Yield occurs when s1 - s3 reaches k (= direct shear strength).
      • This assumes s2 has no effect.
      • Tresca's criterion is more conservative than von Mises's.

    • 3rd. von Mises yield criterion ( called J2-plasticity or J2 flow theory )
      • Yield occurs : 2nd invariant of deviator stress tensor (J2) reaches k^2.
        • von Mises 응력 : 각 지점에서의 비틀림에너지 (max. distortion E) 나타냄.
        • J2 = 1/6 * [(s_xx-s_yy)^2 +...] + s_xy^2 +... = 1/2 * s'_ij * s'_ij
        • material yield strength (k) = direct shear strength
      • 식으로 쓰면 f = J2 - k^2 = 1/2 * (s'_ij)^2 - k^2
        • f = 0 : yield
        • f < 0 : elastic
        • f > 0 : not possible
      • All components of shear stress contribute to yield.

    • 4th. Drucker-Prager yield criterion
      • Pressure-dependent model (for plastic deformation of soils)
      • sqrt(J2) = A + B * I1
        • I1 : 1st invariant of Cauchy stress = 3*mean stress(s_m)
        • J2 : 2nd deviator stress invariant
        • A & B from experiments


    추가 : Hyper-elastic = Green elastic

    s_ij = B_ij + C_ijkl * e_kl

    where, B_ij : components of initial stress tensor = storing of past Energy
    아쉽게도 모든 지반 탄성체는 B_ij 가 당연히 존재함.



    그리하여 


    지반의 불균질 & 비등방성으로 인해, 복잡하지만 탄소성변형을 알아야 함.







    2015년 10월 21일 수요일

    Finite Element Method (유한요소법)



    지반공학을 위한 수치해석, 정말 쉽게 한번 접근해 보자.

    진짜 중요하지만 너무 어려운 유한요소법을 요점정리 해본다.



    대전제: 모두 다 이해하려 하지 마라.



    FEM의 시발점은 간단하게 말하려 해도... shape function 을 기반으로 한 basis function 의 적용과

    Guassian interpolation 에 의한 global system 과 local system 사이의 mapping 이다.

    이게 무슨 소린가 싶겠지만, 아쉽게도 이 부분은 자세히 파헤치지 않는것이 신상에 좋다.

    ㅡㅡㅋ


    Equation of Elasticity: FEM


    [e] = [B] [u(i)]

    [s] = [C] [e] = [C] [B] [u(i)]

    [Ke] = Integral ( [B]' [C] [B] dV )

    [Kg] [u] = [F]



    Process


    1. Initialization
      1. Input values: Geometry, Material, Load
      2. Gaussian interpolation points
      3. Memory allocation
    2. Mesh creation
      1. Global node coordinates
      2. Connectivity matrix: connect from each element to 4 nodes
    3. Stiffness Matrix: Consider each element
      1. [J]: Jacobian matrix
      2. [B]: derivatives of the shape functions matrix
      3. [C]: Constitutive matrix
      4. [Ke] = [B]'[C][B] * wi wj * detJ * t
      5. [Kg]: Assemble local [Ke] to global [Kg]
    4. Boundary Conditions and Load vector
      1. Left
      2. Right
      3. Bottom
      4. [F]: Load vector
    5. Solution
      1. Matrix inversion: [u] = Inv[Kg] * [F]
      2. Post-process for Stress: [s] = [C][B][u] at interpolation points
      3. Interpolate [s] at global nodes




    추가


    Implicit 와 Explicit 의 차이 설명 글


    http://bbalog.tistory.com/12




    추후 다시 업데이트 하겠지만, 우선 간단한 예제를 아래 남긴다.

    Example - Circular Footing on Soil media (axisymmetry)