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
