33 #include "libmesh/parsed_function.h"
34 #include "libmesh/zero_function.h"
59 std::string base_function =
63 if (base_function ==
"0")
64 libmesh_error_msg(
"Error! Zero AveragedFan specified!" <<
67 if (base_function ==
"0")
68 this->base_velocity_function.reset
69 (
new libMesh::ZeroFunction<libMesh::Number>());
71 this->base_velocity_function.reset
72 (
new libMesh::ParsedFunction<libMesh::Number>(base_function));
74 std::string vertical_function =
78 if (vertical_function ==
"0")
79 libmesh_error_msg(
"Warning! Zero LocalVertical specified!" <<
82 this->local_vertical_function.reset
83 (
new libMesh::ParsedFunction<libMesh::Number>(vertical_function));
85 std::string lift_function_string =
89 if (lift_function_string ==
"0")
90 std::cout <<
"Warning! Zero lift function specified!" << std::endl;
92 this->lift_function.reset
93 (
new libMesh::ParsedFunction<libMesh::Number>(lift_function_string));
95 std::string drag_function_string =
99 if (drag_function_string ==
"0")
100 std::cout <<
"Warning! Zero drag function specified!" << std::endl;
102 this->drag_function.reset
103 (
new libMesh::ParsedFunction<libMesh::Number>(drag_function_string));
105 std::string chord_function_string =
109 if (chord_function_string ==
"0")
110 libmesh_error_msg(
"Warning! Zero chord function specified!" <<
113 this->chord_function.reset
114 (
new libMesh::ParsedFunction<libMesh::Number>(chord_function_string));
116 std::string area_function_string =
120 if (area_function_string ==
"0")
121 libmesh_error_msg(
"Warning! Zero area_swept_function specified!" <<
124 this->area_swept_function.reset
125 (
new libMesh::ParsedFunction<libMesh::Number>(area_function_string));
127 std::string aoa_function_string =
129 std::string(
"00000"));
131 if (aoa_function_string ==
"00000")
132 libmesh_error_msg(
"Warning! No angle-of-attack specified!" <<
135 this->aoa_function.reset
136 (
new libMesh::ParsedFunction<libMesh::Number>(aoa_function_string));
141 (
const libMesh::Point& point,
142 const libMesh::Real time,
143 const libMesh::NumberVectorValue& U,
144 libMesh::NumberVectorValue& F,
145 libMesh::NumberTensorValue *dFdU)
148 libmesh_assert(base_velocity_function.get());
150 libMesh::DenseVector<libMesh::Number> output_vec(3);
152 (*base_velocity_function)(point, time,
155 const libMesh::NumberVectorValue U_B(output_vec(0),
159 const libMesh::Number U_B_size = U_B.size();
166 const libMesh::NumberVectorValue N_B =
167 libMesh::NumberVectorValue(U_B/U_B_size);
169 (*local_vertical_function)(point, time,
173 const libMesh::NumberVectorValue N_V(output_vec(0),
179 const libMesh::NumberVectorValue N_R = N_B.cross(N_V);
182 const libMesh::NumberVectorValue U_P = U - (U*N_R)*N_R - U_B;
184 const libMesh::Number U_P_size = U_P.size();
193 const libMesh::NumberVectorValue N_drag =
194 libMesh::NumberVectorValue(-U_P/U_P_size);
197 const libMesh::NumberVectorValue N_lift = N_drag.cross(N_R);
200 const libMesh::Number u_fwd = -(U_P * N_B);
203 const libMesh::Number u_up = U_P * N_V;
207 libmesh_assert (u_up || u_fwd);
210 const libMesh::Number part_angle = std::atan2(u_up, u_fwd);
213 const libMesh::Number angle = part_angle +
214 (*aoa_function)(point, time);
216 const libMesh::Number C_lift = (*lift_function)(point, angle);
217 const libMesh::Number C_drag = (*drag_function)(point, angle);
219 const libMesh::Number chord = (*chord_function)(point, time);
220 const libMesh::Number area = (*area_swept_function)(point, time);
222 const libMesh::Number v_sq = U_P*U_P;
224 const libMesh::Number LDfactor = 0.5 * this->_rho * v_sq * chord / area;
225 const libMesh::Number lift = C_lift * LDfactor;
226 const libMesh::Number drag = C_drag * LDfactor;
229 F = lift * N_lift + drag * N_drag;
235 const libMesh::NumberVectorValue LDderivfactor =
236 (N_lift*C_lift+N_drag*C_drag) *
237 this->_rho * chord / area;
239 for (
unsigned int i=0; i != 3; ++i)
240 for (
unsigned int j=0; j != 3; ++j)
241 (*dFdU)(i,j) = LDderivfactor(i) * U_P(j);
const PhysicsName incompressible_navier_stokes
Physics class for Incompressible Navier-Stokes.
bool compute_force(const libMesh::Point &point, const libMesh::Real time, const libMesh::NumberVectorValue &U, libMesh::NumberVectorValue &F, libMesh::NumberTensorValue *dFdU=NULL)
INSTANTIATE_INC_NS_SUBCLASS(AveragedFanBase)
virtual void read_input_options(const GetPot &input)
Read options from GetPot input file.
const PhysicsName averaged_fan