 Contents of /trunk/doc/user/wave.tex

Revision 4821 - (show annotations)
Tue Apr 1 04:58:33 2014 UTC (5 years, 2 months ago) by sshaw
File MIME type: application/x-tex
File size: 17177 byte(s)
moved SolverOptions to c++, split into SolverOptions for the options and SolverBuddy as the state as a precursor to per-pde solving... does break some use cases (e.g. pde.getSolverOptions().DIRECT will now fail, new value access is with SolverOptions.DIRECT), examples and documentation updated to match

 1 \section{Wave Propagation} 2 \label{WAVE CHAP} 3 4 In this next example we want to calculate the displacement field $u_{i}$ for any time $t>0$ by solving the wave equation:\index{wave equation} 5 \begin{eqnarray}\label{WAVE general problem} 6 \rho u_{i,tt} - \sigma_{ij,j}=0 7 \end{eqnarray} 8 in a three dimensional block of length $L$ in $x_{0}$ 9 and $x_{1}$ direction and height $H$ in $x_{2}$ direction. 10 $\rho$ is the known density which may be a function of its location. 11 $\sigma_{ij}$ is the stress field\index{stress} which in case of an 12 isotropic, linear elastic material is given by 13 \begin{eqnarray}\label{WAVE stress} 14 \sigma_{ij} = \lambda u_{k,k} \delta_{ij} + \mu (u_{i,j} + u_{j,i}) 15 \end{eqnarray} 16 where $\lambda$ and $\mu$ are the Lame coefficients\index{Lame coefficients} 17 and $\delta_{ij}$ denotes the Kronecker symbol\index{Kronecker symbol}. 18 On the boundary the normal stress is given by 19 \begin{eqnarray}\label{WAVE natural} 20 \sigma_{ij}n_{j}=0 21 \end{eqnarray} 22 for all time $t>0$. 23 24 \begin{figure}[t!] 25 \centerline{\includegraphics[width=14cm]{waveu}} 26 \caption{Input Displacement at Source Point ($\alpha=0.7$, $t_{0}=3$, $U_{0}=1$).} 27 \label{WAVE FIG 1.2} 28 \end{figure} 29 30 \begin{figure} 31 \centerline{\includegraphics[width=14cm]{wavea}} 32 \caption{Input Acceleration at Source Point ($\alpha=0.7$, $t_{0}=3$, $U_{0}=1$).} 33 \label{WAVE FIG 1.1} 34 \end{figure} 35 36 Here we are modelling a point source at the point $x_C$ in the 37 $x_{0}$-direction which is a negative pulse of amplitude 38 $U_{0}$ followed by the same positive pulse. 39 In mathematical terms we use 40 \begin{eqnarray} \label{WAVE source} 41 u_{0}(x_C,t)= U_{0} \; \sqrt{2} \frac{(t-t_{0})}{\alpha} e^{\frac{1}{2}-\frac{(t-t_{0})^2}{\alpha^2}} \ 42 \end{eqnarray} 43 for all $t\ge0$ where $\alpha$ is the width of the pulse and $t_{0}$ 44 is the time when the pulse changes from negative to positive. 45 In the simulations we will choose $\alpha=0.3$ and $t_{0}=2$ (see 46 \fig{WAVE FIG 1.2}) and apply the source as a constraint in a sphere of small 47 radius around the point $x_C$. 48 49 We use an explicit time integration scheme to calculate the displacement field 50 $u$ at certain time marks $t^{(n)}$, where $t^{(n)}=t^{(n-1)}+h$ with time 51 step size $h>0$. 52 In the following the upper index ${(n)}$ refers to values at time $t^{(n)}$. 53 We use the Verlet scheme\index{Verlet scheme} with constant time step size $h$ 54 which is defined by 55 \begin{eqnarray} \label{WAVE dyn 2} 56 u^{(n)}=2u^{(n-1)}-u^{(n-2)} + h^2 a^{(n)} \\ 57 \end{eqnarray} 58 for all $n=2,3,\ldots$. It is designed to solve a system of equations of the form 59 \begin{eqnarray} \label{WAVE dyn 2b} 60 u_{,tt}=G(u) 61 \end{eqnarray} 62 where one sets $a^{(n)}=G(u^{(n-1)})$. 63 64 In our case $a^{(n)}$ is given by 65 \begin{eqnarray}\label{WAVE dyn 3} 66 \rho a^{(n)}_{i}=\sigma^{(n-1)}_{ij,j} 67 \end{eqnarray} 68 and boundary conditions 69 \begin{eqnarray} \label{WAVE natural at n} 70 \sigma^{(n-1)}_{ij}n_{j}=0 71 \end{eqnarray} 72 derived from \eqn{WAVE natural} where 73 \begin{eqnarray} \label{WAVE dyn 3a} 74 \sigma_{ij}^{(n-1)} & = & \lambda u^{(n-1)}_{k,k} \delta_{ij} + \mu ( u^{(n-1)}_{i,j} + u^{(n-1)}_{j,i}). 75 \end{eqnarray} 76 We also need to apply the constraint 77 \begin{eqnarray} \label{WAVE dyn 3aa} 78 a^{(n)}_{0}(x_C,t)= U_{0} 79 \frac{\sqrt(2.)}{\alpha^2} (4\frac{(t-t_{0})^3}{\alpha^3}-6\frac{t-t_{0}}{\alpha})e^{\frac{1}{2}-\frac{(t-t_{0})^2}{\alpha^2}} 80 \end{eqnarray} 81 which is derived from equation~\ref{WAVE source} by calculating the second 82 order time derivative (see \fig{WAVE FIG 1.1}). 83 Now we have converted our problem for displacement, $u^{(n)}$, into a problem 84 for acceleration, $a^{(n)}$, which depends on the solution at the previous two 85 time steps $u^{(n-1)}$ and $u^{(n-2)}$. 86 87 In each time step we have to solve this problem to get the acceleration $a^{(n)}$, 88 and we will use the \LinearPDE class of the \linearPDEs package to do so. 89 The general form of the PDE defined through the \LinearPDE class is discussed 90 in \Sec{SEC LinearPDE}. The form which is relevant here is 91 \begin{equation}\label{WAVE dyn 100} 92 D_{ij} a^{(n)}_{j} = - X_{ij,j}\; . 93 \end{equation} 94 The natural boundary condition 95 \begin{equation}\label{WAVE dyn 101} 96 n_{j}X_{ij} =0 97 \end{equation} 98 is used. 99 With $u=a^{(n)}$ we can identify the values to be assigned to $D$ and $X$: 100 \begin{equation}\label{WAVE dyn 6} 101 \begin{array}{ll} 102 D_{ij}=\rho \delta_{ij}& 103 X_{ij}=-\sigma^{(n-1)}_{ij}\\ 104 \end{array} 105 \end{equation} 106 Moreover we need to define the location $r$ where the constraint~\ref{WAVE dyn 3aa} is applied. 107 We will apply the constraint on a small sphere of radius $R$ around 108 $x_C$ (we will use 3\% of the width of the domain): 109 \begin{equation}\label{WAVE dyn 6 1} 110 q_{i}(x) = 111 \left\{ 112 \begin{array}{l l} 113 1 & \quad \text{where $\|x-x_c\|\le R$}\\ 114 0 & \quad \text{otherwise} 115 \end{array} 116 \right. 117 \end{equation} 118 The following script defines the function \function{wavePropagation} which 119 implements the Verlet scheme to solve our wave propagation problem. 120 The argument \var{domain} which is a \Domain class object defines the domain 121 of the problem. 122 \var{h} and \var{tend} are the time step size and the end time of the simulation. 123 \var{lam}, \var{mu} and \var{rho} are material properties. 124 \begin{python} 125 def wavePropagation(domain,h,tend,lam,mu,rho, x_c, src_radius, U0): 126 # lists to collect displacement at point source which is returned 127 # to the caller 128 ts, u_pc0, u_pc1, u_pc2 = [], [], [], [] 129 130 x=domain.getX() 131 # ... open new PDE ... 132 mypde=LinearPDE(domain) 133 mypde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING) 134 kronecker=identity(mypde.getDim()) 135 dunit=numpy.array([1., 0., 0.]) # defines direction of point source 136 mypde.setValue(D=kronecker*rho, q=whereNegative(length(x-xc)-src_radius)*dunit) 137 # ... set initial values .... 138 n=0 139 # for first two time steps 140 u=Vector(0., Solution(domain)) 141 u_last=Vector(0., Solution(domain)) 142 t=0 143 # define the location of the point source 144 L=Locator(domain, xc) 145 # find potential at point source 146 u_pc=L.getValue(u) 147 print("u at point charge=",u_pc) 148 ts.append(t) 149 u_pc0.append(u_pc) 150 u_pc1.append(u_pc) 151 u_pc2.append(u_pc) 152 153 while t

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision