@@ -161,8 +161,8 @@ namespace SPH
161161 size_t index_particle_j = neighboring_particle.j_ ;
162162 BaseParticleData &base_particle_data_j = particles_->base_particle_data_ [index_particle_j];
163163
164- div_correction +=
165- dot (neighboring_particle.gradW_ij_ , neighboring_particle.r_ij_ )* base_particle_data_j.Vol_ ;
164+ div_correction += neighboring_particle. dW_ij_ * neighboring_particle. r_ij_ *
165+ dot (- neighboring_particle.e_ij_ , neighboring_particle.e_ij_ ) * base_particle_data_j.Vol_ ;
166166 }
167167
168168 fluid_data_i.div_correction_ = div_correction;
@@ -183,8 +183,8 @@ namespace SPH
183183 SolidParticleData &solid_data_j
184184 = (*interacting_particles_[interacting_body_index]).solid_body_data_ [index_particle_j];
185185
186- div_correction +=
187- dot (neighboring_particle.gradW_ij_ , neighboring_particle.r_ij_ )*base_particle_data_j.Vol_ ;
186+ div_correction += neighboring_particle. dW_ij_ * neighboring_particle. r_ij_ *
187+ dot (- neighboring_particle.e_ij_ , neighboring_particle.e_ij_ )*base_particle_data_j.Vol_ ;
188188 }
189189
190190 fluid_data_i.div_correction_ += div_correction;
@@ -214,7 +214,7 @@ namespace SPH
214214
215215 // viscous force
216216 Vecd vel_detivative = (base_particle_data_i.vel_n_ - base_particle_data_j.vel_n_ )
217- / (neighboring_particle.r_ij_ . norm () + 0.01 * smoothing_length_);
217+ / (neighboring_particle.r_ij_ + 0.01 * smoothing_length_);
218218 acceleration += 2.0 *mu_*vel_detivative*base_particle_data_j.Vol_ / fluid_data_i.rho_n_
219219 *neighboring_particle.dW_ij_ ;
220220 }
@@ -241,9 +241,9 @@ namespace SPH
241241
242242 // viscous force with a simple wall model for high-Reynolds number flow
243243 Vecd vel_detivative = 2.0 *(base_particle_data_i.vel_n_ - solid_data_j.vel_ave_ )
244- / (neighboring_particle.r_ij_ . norm () + 0.01 * smoothing_length_);
244+ / (neighboring_particle.r_ij_ + 0.01 * smoothing_length_);
245245 Real vel_difference = 0.03 *(base_particle_data_i.vel_n_ - solid_data_j.vel_ave_ ).norm ()
246- *neighboring_particle.r_ij_ . norm () ;
246+ * neighboring_particle.r_ij_ ;
247247 acceleration += 2.0 *SMAX (mu_, fluid_data_i.rho_n_ *vel_difference)
248248 *vel_detivative*base_particle_data_j.Vol_ / fluid_data_i.rho_n_
249249 *neighboring_particle.dW_ij_ ;
@@ -267,10 +267,10 @@ namespace SPH
267267 FluidParticleData &fluid_data_j = particles_->fluid_particle_data_ [index_particle_j];
268268
269269 // exra stress
270- acceleration += 0.5 *dt*((fluid_data_i.rho_n_ * base_particle_data_i.vel_n_
271- *dot (fluid_data_i.dvel_dt_trans_ , neighboring_particle.gradW_ij_ ))
270+ acceleration += 0.5 *dt*((fluid_data_i.rho_n_ * base_particle_data_i.vel_n_
271+ * neighboring_particle. dW_ij_ * dot (fluid_data_i.dvel_dt_trans_ , - neighboring_particle.e_ij_ ))
272272 + (fluid_data_j.rho_n_ *base_particle_data_j.vel_n_
273- * dot (fluid_data_j.dvel_dt_trans_ , neighboring_particle.gradW_ij_ )))
273+ * neighboring_particle. dW_ij_ * dot (fluid_data_j.dvel_dt_trans_ , - neighboring_particle.e_ij_ )))
274274 *base_particle_data_j.Vol_ / fluid_data_i.rho_n_ ;
275275 }
276276
@@ -296,7 +296,7 @@ namespace SPH
296296
297297 // exra stress
298298 acceleration += 0.5 *dt*(fluid_data_i.rho_n_ *base_particle_data_i.vel_n_
299- *dot (fluid_data_i.dvel_dt_trans_ , neighboring_particle.gradW_ij_ ))
299+ * neighboring_particle. dW_ij_ * dot (fluid_data_i.dvel_dt_trans_ , - neighboring_particle.e_ij_ ))
300300 *base_particle_data_j.Vol_ / fluid_data_i.rho_n_ ;
301301 }
302302
@@ -325,8 +325,8 @@ namespace SPH
325325 FluidParticleData &fluid_data_j = particles_->fluid_particle_data_ [index_particle_j];
326326
327327 // acceleration for transport velocity
328- acceleration_trans -= 2.0 * p_background_*base_particle_data_j.Vol_ / fluid_data_i.rho_n_
329- *neighboring_particle.gradW_ij_ ;
328+ acceleration_trans -= 2.0 * p_background_*base_particle_data_j.Vol_ / fluid_data_i.rho_n_
329+ * neighboring_particle.dW_ij_ * (-neighboring_particle. e_ij_ ) ;
330330 }
331331 fluid_data_i.dvel_dt_trans_ = acceleration_trans;
332332 base_particle_data_i.pos_n_ += acceleration_trans * dt*dt*0.5 ;
@@ -348,7 +348,7 @@ namespace SPH
348348
349349 // acceleration for transport velocity
350350 acceleration_trans -= 2.0 *p_background_*base_particle_data_j.Vol_ / fluid_data_i.rho_n_
351- *neighboring_particle.gradW_ij_ ;
351+ * neighboring_particle.dW_ij_ * (-neighboring_particle. e_ij_ ) ;
352352 }
353353
354354 fluid_data_i.dvel_dt_trans_ += acceleration_trans;
@@ -455,7 +455,7 @@ namespace SPH
455455 Real p_star = material_->RiemannSolverForPressure (fluid_data_i.rho_n_ , fluid_data_j.rho_n_ ,
456456 fluid_data_i.p_ , fluid_data_j.p_ , ul, ur);
457457 acceleration -= 2.0 *p_star*base_particle_data_j.Vol_ / fluid_data_i.rho_n_
458- *neighboring_particle.gradW_ij_ ;
458+ * neighboring_particle.dW_ij_ * (-neighboring_particle. e_ij_ ) ;
459459 }
460460
461461 base_particle_data_i.dvel_dt_ = acceleration;
@@ -481,16 +481,17 @@ namespace SPH
481481 Real face_wall_external_acceleration
482482 = dot ((external_force_->InducedAcceleration (base_particle_data_i.pos_n_ , base_particle_data_i.vel_n_ , dt)
483483 - solid_data_j.dvel_dt_ave_ ), -solid_data_j.n_ );
484- Real p_in_wall = fluid_data_i.p_ + fluid_data_i.rho_n_ *dot (neighboring_particle.r_ij_ , -solid_data_j.n_ )
485- *SMAX (0.0 , face_wall_external_acceleration);
484+ Real p_in_wall = fluid_data_i.p_ + fluid_data_i.rho_n_ * neighboring_particle.r_ij_ *
485+ dot (neighboring_particle.e_ij_ , -solid_data_j.n_ )
486+ *SMAX (0.0 , face_wall_external_acceleration);
486487 Real rho_in_wall = material_->ReinitializeRho (p_in_wall);
487488
488489 Real p_star = (fluid_data_i.p_ *rho_in_wall + p_in_wall * fluid_data_i.rho_n_ )
489490 / (fluid_data_i.rho_n_ + rho_in_wall);
490491
491492 // pressure force
492493 acceleration -= 2.0 *p_star*base_particle_data_j.Vol_ / fluid_data_i.rho_n_
493- *neighboring_particle.gradW_ij_ ;
494+ * neighboring_particle.dW_ij_ * (-neighboring_particle. e_ij_ ) ;
494495 }
495496
496497 base_particle_data_i.dvel_dt_ += acceleration;
@@ -552,7 +553,8 @@ namespace SPH
552553 // seems not using rimann problem solution is better
553554 Vecd vel_in_wall = 2.0 *solid_data_j.vel_ave_ - base_particle_data_i.vel_n_ ;
554555 density_change_rate += fluid_data_i.rho_n_ *base_particle_data_j.Vol_
555- *dot (base_particle_data_i.vel_n_ - vel_in_wall, neighboring_particle.gradW_ij_ );
556+ * dot (base_particle_data_i.vel_n_ - vel_in_wall, -neighboring_particle.e_ij_ )
557+ * neighboring_particle.dW_ij_ ;
556558 }
557559
558560 fluid_data_i.drho_dt_ += fluid_data_i.div_correction_ *density_change_rate;
@@ -584,7 +586,7 @@ namespace SPH
584586
585587 // pressure force
586588 acceleration -= 2.0 *p_star*base_particle_data_j.Vol_ / fluid_data_i.rho_n_
587- *neighboring_particle.gradW_ij_ ;
589+ * neighboring_particle.dW_ij_ * (-neighboring_particle. e_ij_ ) ;
588590 }
589591
590592 base_particle_data_i.dvel_dt_ = acceleration;
@@ -610,7 +612,8 @@ namespace SPH
610612 Real face_wall_external_acceleration
611613 = dot ((external_force_->InducedAcceleration (base_particle_data_i.pos_n_ , base_particle_data_i.vel_n_ , dt)
612614 - solid_data_j.dvel_dt_ave_ ), -solid_data_j.n_ );
613- Real p_in_wall = fluid_data_i.p_ + fluid_data_i.rho_n_ *dot (neighboring_particle.r_ij_ , -solid_data_j.n_ )
615+ Real p_in_wall = fluid_data_i.p_ + fluid_data_i.rho_n_ * neighboring_particle.r_ij_ *
616+ dot (neighboring_particle.e_ij_ , -solid_data_j.n_ )
614617 *SMAX (0.0 , face_wall_external_acceleration);
615618 Real rho_in_wall = material_->ReinitializeRho (p_in_wall);
616619
@@ -619,7 +622,7 @@ namespace SPH
619622
620623 // pressure force
621624 acceleration -= 2.0 *p_star*base_particle_data_j.Vol_ / fluid_data_i.rho_n_
622- *neighboring_particle.gradW_ij_ ;
625+ * neighboring_particle.dW_ij_ * (-neighboring_particle. e_ij_ ) ;
623626 }
624627
625628 base_particle_data_i.dvel_dt_ += acceleration;
@@ -693,8 +696,9 @@ namespace SPH
693696 / (fluid_data_i.rho_n_ + fluid_data_j.rho_n_ );
694697
695698 // pressure and elastic force
696- acceleration -= (2.0 *p_star*neighboring_particle.gradW_ij_
697- + (non_newtonian_fluid_data_i.tau_ + non_newtonian_fluid_data_j.tau_ )*neighboring_particle.gradW_ij_ )
699+ acceleration -= (2.0 * p_star * neighboring_particle.dW_ij_ * (-neighboring_particle.e_ij_ )
700+ + (non_newtonian_fluid_data_i.tau_ + non_newtonian_fluid_data_j.tau_ )
701+ * neighboring_particle.dW_ij_ * (-neighboring_particle.e_ij_ ))
698702 *base_particle_data_j.Vol_ / fluid_data_i.rho_n_ ;
699703 }
700704
@@ -722,16 +726,17 @@ namespace SPH
722726 Real face_wall_external_acceleration
723727 = dot ((external_force_->InducedAcceleration (base_particle_data_i.pos_n_ , base_particle_data_i.vel_n_ , dt)
724728 - solid_data_j.dvel_dt_ave_ ), -solid_data_j.n_ );
725- Real p_in_wall = fluid_data_i.p_ + fluid_data_i.rho_n_ *dot (neighboring_particle.r_ij_ , -solid_data_j.n_ )
729+ Real p_in_wall = fluid_data_i.p_ + fluid_data_i.rho_n_ * neighboring_particle.r_ij_ *
730+ dot (neighboring_particle.e_ij_ , -solid_data_j.n_ )
726731 *SMAX (0.0 , face_wall_external_acceleration);
727732 Real rho_in_wall = material_->ReinitializeRho (p_in_wall);
728733
729734 Real p_star = (fluid_data_i.p_ *rho_in_wall + p_in_wall * fluid_data_i.rho_n_ )
730735 / (fluid_data_i.rho_n_ + rho_in_wall);
731736
732737 // stress boundary condition
733- acceleration -= (2.0 * p_star* neighboring_particle.gradW_ij_
734- + 2.0 * non_newtonian_fluid_data_i.tau_ * neighboring_particle.gradW_ij_ )
738+ acceleration -= (2.0 * p_star * neighboring_particle.dW_ij_ * (-neighboring_particle. e_ij_ )
739+ + 2.0 * non_newtonian_fluid_data_i.tau_ * neighboring_particle.dW_ij_ * (-neighboring_particle. e_ij_ ) )
735740 *base_particle_data_j.Vol_ / fluid_data_i.rho_n_ ;
736741 }
737742
@@ -773,7 +778,7 @@ namespace SPH
773778 *(u_star - ul)*neighboring_particle.dW_ij_ ;
774779 // transport velocity
775780 Matd velocity_gradient = SimTK::outer ((fluid_data_i.vel_trans_ - fluid_data_j.vel_trans_ ),
776- neighboring_particle.gradW_ij_ )*base_particle_data_j.Vol_ ;
781+ neighboring_particle.dW_ij_ * (-neighboring_particle. e_ij_ ) )*base_particle_data_j.Vol_ ;
777782 stress_rate -= ~velocity_gradient*non_newtonian_fluid_data_i.tau_
778783 + non_newtonian_fluid_data_i.tau_ *velocity_gradient + non_newtonian_fluid_data_i.tau_ /lambda_
779784 + (~velocity_gradient + velocity_gradient)*mu_p_/lambda_;
@@ -803,10 +808,10 @@ namespace SPH
803808 = (*interacting_particles_[interacting_body_index]).solid_body_data_ [index_particle_j];
804809
805810 density_change_rate += 2.0 *fluid_data_i.rho_n_
806- *dot (fluid_data_i.vel_trans_ - solid_data_j.vel_ave_ , neighboring_particle. gradW_ij_ )
807- * base_particle_data_j.Vol_ ;
811+ *dot (fluid_data_i.vel_trans_ - solid_data_j.vel_ave_ ,
812+ neighboring_particle. dW_ij_ * (-neighboring_particle. e_ij_ )) * base_particle_data_j.Vol_ ;
808813 Matd velocity_gradient = SimTK::outer ((fluid_data_i.vel_trans_ - solid_data_j.vel_ave_ ),
809- neighboring_particle.gradW_ij_ )* base_particle_data_j.Vol_ *2.0 ;
814+ neighboring_particle.dW_ij_ * (-neighboring_particle. e_ij_ )) * base_particle_data_j.Vol_ *2.0 ;
810815 stress_rate -= ~velocity_gradient*non_newtonian_fluid_data_i.tau_
811816 + non_newtonian_fluid_data_i.tau_ *velocity_gradient + non_newtonian_fluid_data_i.tau_ / lambda_
812817 + (~velocity_gradient + velocity_gradient)*mu_p_ / lambda_;
0 commit comments