Changes between Version 5 and Version 6 of u/johannjc/injection


Ignore:
Timestamp:
06/07/18 12:42:40 (7 years ago)
Author:
Jonathan
Comment:

Legend:

Unmodified
Added
Removed
Modified
  • u/johannjc/injection

    v5 v6  
    88First, 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
    99
    10 $\delta m = M K(x)$
     10$m_i = M K( |\vec{x}_i-\vec{X} |)$
    1111
    1212where our normalized kernel
    1313
    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$.
    1515
    1616Here $r$ is the distance from a sink particle.
     
    2222Note 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.
    2323
    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
     24We 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
     32If we consider each zone inside the kernel, we can write the constraints as the following matrix equation
    2533
    2634{{{
     
    2836 \[
    2937 \left[\begin{array}{rrrrr}
    30   x_1 & x_2 & x_3 & ... & x_n  \\
    31    y_1 & y_2 & y_3 & ... & y_n  \\
    3238   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}  \\
    3341   \end{array} \right]
    3442 \left[\begin{array}{r}
    35   \delta m_1  \\
    36   \delta m_2  \\
    37   \delta m_3  \\
     43  m_1  \\
     44  m_2  \\
     45  m_3  \\
    3846  ...  \\
    39   \delta m_n  \\
     47  m_n  \\
    4048   \end{array} \right]=
    4149 \left[\begin{array}{r}
    42   MX_0  \\
    43   MY_0  \\
    4450  M  \\
     51  MX_1  \\
     52  MX_2  \\
     53
    4554   \end{array} \right]
    4655\]
     
    5261[[Image(Screen Shot 2017-12-15 at 3.07.00 PM.png​)]]
    5362
    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.
     63While 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.
    5564
    5665We 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)
     
    6271
    6372
    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$
     73Now 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$.
    6574
     75However, 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 
    6677{{{
    6778#!latex
     
    7182  0 & ... & 0 & m_1 & ... & m_n & 0 & ... & 0 \\
    7283  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  \\
    7387   0 & ... & 0  & m_1z_1 & ... & m_nz_n &  -m_1y_1 & ... & -m_ny_n\\
    7488  -m_1z_1 & ... & -m_nz_n &  0 & ... & 0 & m_1x_1 & ... & m_nx_n \\
    7589m_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  \\
    7990   \end{array} \right]
    8091\]
     
    98109  0  \\
    99110  0  \\
     111  P_x \\
     112  P_y \\
     113  P_z \\
    100114  J_x  \\
    101115  J_y  \\
    102116  J_z  \\
    103   P_x \\
    104   P_y \\
    105   P_z \\
    106117   \end{array} \right]
    107118\]
     
    110121where $\hat{x}_i = \frac{x_i}{|x_i|}$
    111122
     123
     124
     125== Taking into account particle motion and orientation ==
     126If we are explicit about our constraints, we have
     127
     128
     129Change in mass of particle
     130
     131$\displaystyle \sum_i m_i = M$
     132
     133
     134Change in mass moment of particle
     135
     136$\displaystyle\sum_i m_i x_{i,j} = M X_j$
     137
     138Change in momentum of particle (due to loss of mass)
     139
     140$\displaystyle\sum_i m_i v_{i,j} = M V_j$
     141
     142Change 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
     146Amount 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
     150which 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
     155An 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