zoukankan      html  css  js  c++  java
  • 1D RKDG to shallow water equations

    RKDG to shallow water equations

    1.Governing Equations

    [frac{partial U}{partial t} + frac{partial F}{partial x} = 0 ]

    [U = egin{bmatrix} h cr q end{bmatrix} quad F = egin{bmatrix} q cr gh^2/2 + q^2/h end{bmatrix}]

    2.Discrete with DGM

    [egin{equation} U_h = sum{l_j U_j} quad F_h(U) = sum{l_j F(U_j)} end{equation} ]

    [egin{equation}int_{Omega} l_i l_j frac{partial U_j}{partial t} dx+ int_{Omega} l_i frac{partial l_j}{partial x} F(U_j) dx= 0 end{equation}]

    [egin{equation} int_{Omega} l_i l_j frac{partial U_j}{partial t} dx + int_{Omega} l_i frac{partial l_j}{partial x} F(U_j) dx+ oint_{partial Omega} l_i l_j (F^* - F)cdot vec{n} ds = 0 end{equation}]

    [egin{equation} JM frac{partial U}{partial t} + JMD_x F(U) + J_E M_E (F^* - F)cdot vec{n} = 0 end{equation} ]

    ODE:

    [egin{equation} frac{partial U}{partial t} = -frac{partial r}{partial x}D_r F(U) + frac{J_E}{J}M^{-1} M_E (F^* - F)cdot vec{n}=L(U(t)) end{equation} ]

    [egin{equation} rhs = -frac{partial r}{partial x}D_r F(U) + frac{J_E}{J}M^{-1} M_E (F - F^*)cdot vec{n}end{equation} ]

    It is important to point out that at dry cells no flux is flow inside the elemnt. Therefor, for dry cells

    [egin{equation} rhs = frac{J_E}{J}M^{-1} M_E (F - F^*)cdot vec{n}end{equation} ]

    3.Numerical Flux

    3.1.HLL flux function

    Formulations are given as

    [F^{HLL} = left{ egin{matrix} F^- cr frac{S_R F^- - S_L F^+ + S_L S_R(U^+ - U^-)}{S_R S_L} cr F^+ end{matrix} ight. egin{matrix} S_L geq 0 cr S_L < 0 < S_R cr S_R leq 0 end{matrix}]

    Wave Speed is suggested by Fraccarollo and Toro (1995)

    [S_L = min(u^- - sqrt{gh^-}, u^* - c^*) ]

    [S_R = min(u^+ + sqrt{gh^+}, u^* + c^*) ]

    (u^*) and (c^*) is defined by

    [u^* = frac{1}{2}(u^- + u^+) + sqrt{gh^-} - sqrt{gh^+} ]

    [c^* = frac{1}{2}(sqrt{gh^-} + sqrt{gh^+}) + frac{1}{4}(u^- - u^+) ]

    for wet-dry interface, the wave speed is giving as

    1. left-hand dry bed

    [egin{equation} S_L = u^+ - 2sqrt{g h^+} quad S_R = u^+ + sqrt{g h^+} end{equation}]

    1. right-hand dry bed

    [egin{equation} S_L = u^- - sqrt{g h^-} quad S_R = u^- + 2sqrt{g h^-} end{equation}]

    1. both sides are dry

    [egin{equation} S_L = 0 quad S_R = 0 end{equation}]

    Noticing. 1
    For flux terms, the discharge (q^2) is divided by water depth (h)

    [F = egin{bmatrix} q cr gh^2/2 + q^2/h end{bmatrix} ]

    so a threadhold of water depth (h_{flux}) ( (10^{-3})m ) is add into flux function SWEFlux.m. When (h) is less than (h_{flux}), the (q^2/h) is approximated to 0 as there is no flow at this node.

    Noticing. 2
    When defining the dry beds, another threadhold of water depth (h_{dry}) is used. It is convenient to deine (h_{dry}) equals to (h_{flux}).

    3.2.Rotational invariance

    [T = egin{bmatrix} 1 & 0 cr 0 & n_xend{bmatrix} quad T^{-1} = egin{bmatrix} 1 & 0 cr 0 & n_xend{bmatrix}]

    [mathbf{F} cdot mathbf{n} = mathbf{F} cdot n_x = T^{-1}mathbf{F}(TU) ]

    defining (Q = TU), the numerical flux (hat{mathbf{F}}) can be obtained through the evaluation of numerical flux (mathbf{F}) by

    [hat{mathbf{F}} cdot n = T^{-1}mathbf{F}^{HLL}(Q) ]

    4.Limiter

    Note: discontinuity detector from Krivodonova (2003) is not working

    For better numerical stability, minmod limiter is used in limiting the discharge and elevation.

    Check testing/Limiter1D/doc for more details about minmod limiter.

    5. Positive preserving limiter

    For the thin layer approach, a small depth ( (h_{positive} = 10^{-3} m)) and zeros velocity are prescribed for dry nodes.

    The first step is to define wet elements. After each time step, the whole domain is calculated; If the any depth of nodes in (Omega_i) is greater than (h_{positive}), then the element is defined as wet element, otherwise the water height of all nodes are remain unchanged.

    The second step is to modify wet cells; If the depth of any nodes is less than (h_{positive}), then the flow rate is reset to zero and the new water depth is constructed as

    [egin{equation} mathrm{M}Pi_h h_i(x) = heta_1 left( h_i(x) - ar{h}_i ight) + ar{h}_i end{equation}]

    where

    [egin{equation} heta_1 = min left{ frac{ar{h}_i - xi }{ar{h}_i - h_{min}}, 1 ight}, quad h_{min} = min{ h_i (x_i) } end{equation}]

    It is necessary to fulfill the restriction that the mean depth (ar{h}_i) is greater than (xi), i.e. (10^{-4})m. In the function PositiveOperator, if the mean depth of element is less than (xi), all nodes will add a small depth (xi - ar{h}_i) to re-fulfill the restriction.

    At last, all values of water height at nodes with negative (h_i(x_j) <0) will be modified to zero and the discharge of dry nodes ( (h_i le h_{positive}) ) will be reseted to zero.

    6. Wet/Dry reconstruction

    No special treatment is introduced in the model at the moment.

    5.Numerical Test

    5.1.Wet dam break

    Model Setting value
    channel length 1000m
    dam position 500m
    upstream depth 10m
    downstream depth 2m
    element num 400
    Final Time 20s

    5.2.Dry dam break

    Model Setting value
    channel length 1000m
    dam position 500m
    upstream depth 10m
    downstream depth 0m
    element num 400
    Final Time 20s

    5.3.Parabolic bowl

    Model Setting value
    channel length 2000m
    (h_0) 10m
    (a) 600m
    (B) 5m/s
    (T) 269s

    Exact solution

    [egin{equation} Z(x,t) = frac{-B^2 mathrm{cos}(2wt) - B^2 - 4Bw mathrm{cos}(wt)x}{4g} end{equation}]

    1. (t = T/2)

    2. (t = 3T/4)

    3. (t = T)

  • 相关阅读:
    mysql 中索引的限制
    lvm扩展记录
    转载:权威GIS类数据网站汇总
    转载:文件系统inodes使用率过高问题处理
    转载: k8s--pod的状态为evicted
    转载:k8s更新策略
    转载:Tomcat的JVM内存溢出解决方法
    软件推荐
    U盘安装Centos7 问题记录
    转载:Linux下查找文件
  • 原文地址:https://www.cnblogs.com/li12242/p/5317387.html
Copyright © 2011-2022 走看看