Hello Zhang,
What happens if you just try to compute the viscous
drag (or vice versa) ? If one of the drags gives a reasonable answer and
the other diverges, we can narrow things down.
Thanks.
On Sun, Apr 26, 2015 at 11:06 AM, Zhang <[email protected]> wrote:
> Hi,
>
> I got some problems when writing some functions for normal, axial forces
> and pitching moments of an airfoil in an incompressible laminar flow
> (Re=68.8).
>
> Here I calculate the forces on the wing boundary elements as follows,
> \vec{F}=\sum_{qpface_bound} {\tau*normal -normal*p_{qp}}
> *Jxw_face_bound[qp]
> \vec{M}=\sum_{qpface_bound} R.cross({\tau*normal -normal*p_{qp}})
> *Jxw_face_bound[qp]
> where
> p_{qp}=\sum_{0<=m<n_psi_face}{psi_face[m][qp]*soln(dof_p[m])},
> normal=(fe_vel_face->get_normals())[qp];
> tau=2*S*visc;
> S=0.5*(gradu[m](n)+gradu[n](m)), m,n=0,1,2
> R=fe_vel_face->get_xyz()-center_xyz;
>
> When I run the code, sometimes it gives quite large nonphysical force
> components values like -4.698877e+256 , even at serial run. I guess I must
> have made some improper usage of face quadratures.
> Should I use only the velocity quadrature, or with the pressure one as
> well?
>
> Here is the original code section for this work:
>
> template <unsigned int dimension>
> void BoundaryMoment<dimension>::compute()
> {
> EquationSystems &es = *_equation_system;
>
> const MeshBase& mesh = es.get_mesh();
>
> const unsigned int dim = mesh.mesh_dimension();
>
> TransientNonlinearImplicitSystem& system =
> es.get_system<TransientNonlinearImplicitSystem>(this->system_name);
>
> NumericVector<Number>& soln =*(system.current_local_solution);
> const unsigned int u_var = system.variable_number ("u");
> const unsigned int p_var = system.variable_number ("p");
>
> FEType fe_vel_type = system.variable_type(u_var);
> FEType fe_pres_type = system.variable_type(p_var);
>
> const DofMap& dof_map = system.get_dof_map();
>
> AutoPtr<FEBase> fe_vel_face (FEBase::build(dim, fe_vel_type));
> AutoPtr<FEBase> fe_pres_face (FEBase::build(dim, fe_pres_type));
>
> QGauss qface(dim-1, fe_vel_type.default_quadrature_order());
> QGauss qface_pres(dim-1, fe_pres_type.default_quadrature_order());
>
> fe_vel_face->attach_quadrature_rule (&qface);
> fe_pres_face->attach_quadrature_rule (&qface_pres);
>
> const std::vector<std::vector<Real> >& phi_face =
> fe_vel_face->get_phi();
> const std::vector<std::vector<RealGradient> >& dphi_face =
> fe_vel_face->get_dphi();
> const std::vector<Real>& JxW_face = fe_vel_face->get_JxW();
> const std::vector<Point>& qface_normals = fe_vel_face->get_normals();
> const std::vector<Point>& qface_points = fe_vel_face->get_xyz();
>
> const std::vector<std::vector<Real> >& psi_face =
> fe_pres_face->get_phi();
> const std::vector<std::vector<RealGradient> >& dpsi_face =
> fe_pres_face->get_dphi();
> const std::vector<Real>& JxW_face_pres = fe_pres_face->get_JxW();
> const std::vector<Point>& qface_normals_pres =
> fe_pres_face->get_normals();
> const std::vector<Point>& qface_points_pres = fe_pres_face->get_xyz();
>
> std::vector<dof_id_type> dof_indices;
> std::vector<std::vector<dof_id_type> > dof_indices_u;
> dof_indices_u.resize(this->_dim);
> std::vector<dof_id_type> dof_indices_p;
>
> const Real visc = es.parameters.get<Real>("viscosity");
>
> this->_values.zero();
>
> const BoundaryInfo& b_info=mesh.get_boundary_info();
>
> RealTensorValue _tau,_S2;
>
> Number p;
> Number jxw;
>
> Point r; //vector from _origin to any qp point
>
> MeshBase::const_element_iterator el =
> mesh.active_local_elements_begin();
> const MeshBase::const_element_iterator end_el =
> mesh.active_local_elements_end();
>
> int nelem=0;
> for ( ; el != end_el; ++el)
> {
> const Elem* elem = *el;
>
> for(unsigned int d=0;d<this->_dim;d++)
> dof_map.dof_indices (elem, dof_indices_u[d], u_var+d);
> dof_map.dof_indices (elem, dof_indices_p, p_var);
>
> for (unsigned int side=0; side<elem->n_sides(); side++)
> {
>
> bool boundary_found=false;
> for(boundary_id_type& bid:this->_boundary_ids)
> {
> if(b_info.has_boundary_id(elem,side,bid))
> {
> boundary_found=true;
> break;
> }
> }
> if(boundary_found)
> {
> AutoPtr<Elem> _side (elem->build_side(side));
> fe_vel_face->reinit(elem, side);
> fe_pres_face->reinit(elem, side);
> const std::vector<Point>& normals =
> fe_vel_face->get_normals();
>
> const std::vector<Point>&
> qface_points=fe_vel_face->get_xyz(); // physical xyz
>
> for (unsigned int qp=0; qp<qface.n_points(); qp++)
> {
> Gradient gradu[3];
> r=qface_points[qp]-_origin;
> const Point& _normal=normals[qp];
> jxw=JxW_face[qp];
>
> gradu[0].zero();
> gradu[1].zero();
> gradu[2].zero();
> for (unsigned int l=0; l<phi_face.size(); l++)
> {
> for(unsigned int d=0;d<this->_dim;d++)
> gradu[d].add_scaled (dphi_face[l][qp],soln
> (dof_indices_u[d][l]));
> }
>
> _tau.zero();
> _S2.zero();
>
> for(unsigned int m=0; m<this->_dim; m++)
> for(unsigned int n=0; n<this->_dim; n++)
> {
> _S2(m,n)=gradu[m](n);
> }
> _S2+=_S2.transpose();
> _tau=_S2*visc;
>
> p=0.;
> // for (unsigned int l=0; l<n_p_dofs; l++)
> for (unsigned int l=0; l<psi_face.size(); l++)//???
> {
> p += psi_face[l][qp]*soln(dof_indices_p[l]);
> }
>
> this->_values+=r.cross(_tau*_normal-_normal*p)*jxw;
> }// qp loop in qface
> nelem++;
> } // if (elem->on_boundary())
> }// for (unsigned int side=0; side<elem->n_sides(); side++)
> } // el loop
>
> #if 1
> system.comm().sum<Number>(this->_values(0));
> system.comm().sum<Number>(this->_values(1));
> system.comm().sum<Number>(this->_values(2));
> #endif
>
> cout<<"nelem=: "<<nelem<<"; _values: "
> <<this->_values(0)<<","
> <<this->_values(1)<<","
> <<this->_values(2)
> <<endl;
>
> } // void BoundaryMoment::compute()
>
> Thank you for any possible help,
>
> Zhenyu
>
> ------------------------------------------------------------------------------
> One dashboard for servers and applications across Physical-Virtual-Cloud
> Widest out-of-the-box monitoring support with 50+ applications
> Performance metrics, stats and reports that give you Actionable Insights
> Deep dive visibility with transaction tracing using APM Insight.
> http://ad.doubleclick.net/ddm/clk/290420510;117567292;y
> _______________________________________________
> Libmesh-users mailing list
> [email protected]
> https://lists.sourceforge.net/lists/listinfo/libmesh-users
>
--
Vikram Garg
Postdoctoral Associate
Predictive Engineering and Computational Science (PECOS)
The University of Texas at Austin
http://web.mit.edu/vikramvg/www/
http://vikramvgarg.wordpress.com/
http://www.runforindia.org/runners/vikramg
------------------------------------------------------------------------------
One dashboard for servers and applications across Physical-Virtual-Cloud
Widest out-of-the-box monitoring support with 50+ applications
Performance metrics, stats and reports that give you Actionable Insights
Deep dive visibility with transaction tracing using APM Insight.
http://ad.doubleclick.net/ddm/clk/290420510;117567292;y
_______________________________________________
Libmesh-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-users
