Changes between Version 5 and Version 6 of u/johannjc/injection
- Timestamp:
- 06/07/18 12:42:40 (7 years ago)
Legend:
- Unmodified
- Added
- Removed
- Modified
-
u/johannjc/injection
v5 v6 8 8 First, let's consider a simpler example. Let's assume we want to inject a certain amount of mass $M_0$ around a point $\vec{X}$ (in 2D) with a kernel given by 9 9 10 $ \delta m = M K(x)$10 $m_i = M K( |\vec{x}_i-\vec{X} |)$ 11 11 12 12 where our normalized kernel 13 13 14 $K( x) = \frac{\exp(-r^2)}{\pi(1-1/e)}$ for $r < 1$.14 $K(r) = \frac{\exp(-r^2)}{\pi(1-1/e)}$ for $r < 1$. 15 15 16 16 Here $r$ is the distance from a sink particle. … … 22 22 Note the center of mass of the injected material is not .05,.05 and the injected mass is not 1.0. This is because of discretization error. 23 23 24 We would like to find a solution for $\delta m_i$ that is close to this, but subject to the constraints involving the center of mass and the total mass injected. If we consider each zone inside the kernel, we can write the matrix equation 24 We would like to find a solution for $m_i$ that is close to this, but subject to the constraints involving the center of mass and the total mass injected. 25 {{{ 26 #!latex 27 $\sum m_i=M$ 28 29 $\sum m_i \vec{x_i} = M \vec{X}$ 30 }}} 31 32 If we consider each zone inside the kernel, we can write the constraints as the following matrix equation 25 33 26 34 {{{ … … 28 36 \[ 29 37 \left[\begin{array}{rrrrr} 30 x_1 & x_2 & x_3 & ... & x_n \\31 y_1 & y_2 & y_3 & ... & y_n \\32 38 1 & 1 & 1 & ... & 1 \\ 39 x_{1,1} & x_{2,1} & x_{3,1} & ... & x_{n,1} \\ 40 x_{1,2} & x_{2,2} & x_{3,2} & ... & x_{n,2} \\ 33 41 \end{array} \right] 34 42 \left[\begin{array}{r} 35 \deltam_1 \\36 \deltam_2 \\37 \deltam_3 \\43 m_1 \\ 44 m_2 \\ 45 m_3 \\ 38 46 ... \\ 39 \deltam_n \\47 m_n \\ 40 48 \end{array} \right]= 41 49 \left[\begin{array}{r} 42 MX_0 \\43 MY_0 \\44 50 M \\ 51 MX_1 \\ 52 MX_2 \\ 53 45 54 \end{array} \right] 46 55 \] … … 52 61 [[Image(Screen Shot 2017-12-15 at 3.07.00 PM.png)]] 53 62 54 While this solves our constraints - it looks nothing like our desired profile. What we would like to do is to solve $Ax=b$, in a way that does not minimize $||x||$, but that minimizes $||x-d||$ where $d$ is our vector of initial $ \delta_m$ given by our kernel.63 While this solves our constraints - it looks nothing like our desired profile. What we would like to do is to solve $Ax=b$, in a way that does not minimize $||x||$, but that minimizes $||x-d||$ where $d$ is our vector of initial $m_i$ given by our kernel. 55 64 56 65 We can do this using Matlab's lsqlin function - or we can use the prescription given in https://see.stanford.edu/materials/lsoeldsee263/08-min-norm.pdf (page 8-13 through 8-15) … … 62 71 63 72 64 Now for the more complicated case, we have $\delta m_i, v_i$, but the matrix would contain two blocks - so we can solve the problem in two steps. First we can solve for $\delta m_i$, and then use that solution to construct the matrix equation for $v_i$73 Now for the more complicated case, we have terms in involving momentum, and we could simultaneously solve for $m_i$ and $\vec{p}_i$, though it is probably more important to have $v_i$ be close to the desired values. Errors in $m_i$ and $p_i$ could lead to unphysical velocities. We cannot solve for $v_i$ and $m_i$ at the same time, since there are non-linear terms that involve the product of $m_i$ and $v_i$. 65 74 75 However, if we separate out the solution for $m_i$ using the constraints on total mass and mass moment, we can treat those as constants when solving for $v_i$. This then gives us the matrix equation 76 66 77 {{{ 67 78 #!latex … … 71 82 0 & ... & 0 & m_1 & ... & m_n & 0 & ... & 0 \\ 72 83 0 & ... & 0 & 0 & ... & 0 & m_1 & ... & m_n \\ 84 m_1 \hat{x}_1& ... & m_n \hat{x}_n & 0 & ... & 0 & 0 & ... & 0 \\ 85 0 & ... & 0 & m_1\hat{y}_1 & ... & m_n\hat{y}_n & 0 & ... & 0 \\ 86 0 & ... & 0 & 0 & ... & 0 & m_1\hat{z}_1 & ... & m_n\hat{z}_n \\ 73 87 0 & ... & 0 & m_1z_1 & ... & m_nz_n & -m_1y_1 & ... & -m_ny_n\\ 74 88 -m_1z_1 & ... & -m_nz_n & 0 & ... & 0 & m_1x_1 & ... & m_nx_n \\ 75 89 m_1y_1 & ... & m_ny_n & -m_1x_1 & ... & -m_nx_n &0 & ... & 0 \\ 76 m_1 \hat{x}_1& ... & m_n \hat{x}_n & 0 & ... & 0 & 0 & ... & 0 \\77 0 & ... & 0 & m_1\hat{y}_1 & ... & m_n\hat{y}_n & 0 & ... & 0 \\78 0 & ... & 0 & 0 & ... & 0 & m_1\hat{z}_1 & ... & m_n\hat{z}_n \\79 90 \end{array} \right] 80 91 \] … … 98 109 0 \\ 99 110 0 \\ 111 P_x \\ 112 P_y \\ 113 P_z \\ 100 114 J_x \\ 101 115 J_y \\ 102 116 J_z \\ 103 P_x \\104 P_y \\105 P_z \\106 117 \end{array} \right] 107 118 \] … … 110 121 where $\hat{x}_i = \frac{x_i}{|x_i|}$ 111 122 123 124 125 == Taking into account particle motion and orientation == 126 If we are explicit about our constraints, we have 127 128 129 Change in mass of particle 130 131 $\displaystyle \sum_i m_i = M$ 132 133 134 Change in mass moment of particle 135 136 $\displaystyle\sum_i m_i x_{i,j} = M X_j$ 137 138 Change in momentum of particle (due to loss of mass) 139 140 $\displaystyle\sum_i m_i v_{i,j} = M V_j$ 141 142 Change in internal angular momentum of particle 143 144 $\displaystyle\sum_{i,k,l} m_i \epsilon_{j,k,l} x_{i,k} v_{i,l} = J_j$ 145 146 Amount of outward momentum in jet 147 148 $\displaystyle \sum_{i,j} m_i \left ( v_{i,j} - V_j \right ) \frac{x_{i,j}-X_j}{| \vec{x}_{i}-\vec{X} |} = P$ 149 150 which can be reduced to 151 152 $\displaystyle \sum_{i,j} m_i \frac{x_{i,j}-X_j}{| \vec{x}_{i}-\vec{X} |}v_{i,j} = P+\sum_{i,j} m_i V_j \frac{x_{i,j}-X_j}{| \vec{x}_{i}-\vec{X} |}$ 153 154 155 An alternative approach is to try to match the amount of outward momentum in the jet for each direction independently. 156 157 158 $\displaystyle \sum_{i} m_i \left ( v_{i,j} - V_j \right ) \frac{x_{i,j}-X_j}{\left | x_{i,j}-X_j \right |} = P_j$ 159 160 161 $\displaystyle \sum_{i} m_i v_{i,j} \frac{x_{i,j}-X_j}{\left | x_{i,j}-X_j \right |} = P_j+\displaystyle \sum_{i} m_i V_{j} \frac{x_{i,j}-X_j}{\left | x_{i,j}-X_j \right |}$ 162