correction equations fluides

marblesbag/src/fluid_problem.c

@@ -215,12 +215,12 @@ static void fluid_problem_assemble_system(FluidProblem *problem, const double *s
                             {dv[0]-v*dc[0]/c,dv[1]-v*dc[1]/c}};
         local_vector[iphi+0] += jw*(
            2*mu*(dphii[0]*tau[0][0]+dphii[1]*(tau[0][1]+tau[1][0])*0.5)
-          +rho*phii*((u-u_old)/dt+(u*tau[0][0]+v*tau[0][1])/c)
+          +rho*phii*((u-u_old)/dt-u/c *dcdt+(u*tau[0][0]+v*tau[0][1])/c)
           -p*dphii[0]-fu*phii
         );
         local_vector[iphi+3] += jw*(
            2*mu*(dphii[0]*(tau[0][1]+tau[1][0])*0.5+dphii[1]*tau[1][1])
-          +rho*phii*((v-v_old)/dt+(u*tau[1][0]+v*tau[1][1])/c)
+          +rho*phii*((v-v_old)/dt-v/c *dcdt+(u*tau[1][0]+v*tau[1][1])/c)
           -p*dphii[1]-fv*phii
         );
 
@@ -233,14 +233,14 @@ static void fluid_problem_assemble_system(FluidProblem *problem, const double *s
           double dphij[2] = {dphi[jphi][0],dphi[jphi][1]};
           double dtau[2] = {dphij[0]-phij*dc[0]/c,dphij[1]-phij*dc[1]/c};
           local_matrix[(iphi+0)*9+jphi+0] += jw*2*mu*(dphii[0]*dtau[0]+dphii[1]*dtau[1]*0.5)  // U U
-                                            +jw*rho*phii*(phij/dt+(phij*tau[0][0]+u*dtau[0]+v*dtau[1])/c);
+                                            +jw*rho*phii*(phij/dt-phij/c*dcdt+(phij*tau[0][0]+u*dtau[0]+v*dtau[1])/c);
           local_matrix[(iphi+0)*9+jphi+3] += jw*2*mu*dphii[1]*dtau[0]*0.5                     // U V
                                             +jw*rho*phii*phij*tau[0][1]/c;
           local_matrix[(iphi+0)*9+jphi+6] +=-jw*dphii[0]*phij;                                // U P
           local_matrix[(iphi+3)*9+jphi+0] += jw*2*mu*dphii[0]*dtau[1]*0.5                     // V U
                                             +jw*rho*phii*phij*tau[1][0]/c;
           local_matrix[(iphi+3)*9+jphi+3] += jw*2*mu*(dphii[0]*dtau[0]*0.5+dphii[1]*dtau[1])  // V V
-                                            +jw*rho*phii*(phij/dt+(phij*tau[1][1]+u*dtau[0]+v*dtau[1])/c);
+                                            +jw*rho*phii*(phij/dt-phij/c*dcdt+(phij*tau[1][1]+u*dtau[0]+v*dtau[1])/c);
           local_matrix[(iphi+3)*9+jphi+6] +=-jw*dphii[1]*phij;                                // V P
           local_matrix[(iphi+6)*9+jphi+0] += jw*phii*dphij[0];                                // P U
           local_matrix[(iphi+6)*9+jphi+3] += jw*phii*dphij[1];                                // P V