40 base_velocity_function(
""),
41 local_vertical_function(
""),
46 area_swept_function(
""),
55 this->_fan_speed_var = system->add_variable(_fan_speed_var_name,
65 system->time_evolving(this->fan_speed_var());
74 this->set_parameter(base_velocity_function, input,
76 this->zero_vector_function);
78 if (base_velocity_function.expression() == this->zero_vector_function)
79 libmesh_error_msg(
"Error! Zero AveragedTurbine specified!" <<
82 this->set_parameter(local_vertical_function, input,
84 this->zero_vector_function);
86 if (local_vertical_function.expression() == this->zero_vector_function)
87 libmesh_error_msg(
"Error! Zero LocalVertical specified!" <<
90 this->set_parameter(lift_function, input,
94 if (lift_function.expression() ==
"0")
95 std::cout <<
"Warning! Zero lift function specified!" << std::endl;
97 this->set_parameter(drag_function, input,
101 if (drag_function.expression() ==
"0")
102 std::cout <<
"Warning! Zero drag function specified!" << std::endl;
104 this->set_parameter(chord_function, input,
108 if (chord_function.expression() ==
"0")
109 libmesh_error_msg(
"Error! Zero chord function specified!" <<
112 this->set_parameter(area_swept_function, input,
116 if (area_swept_function.expression() ==
"0")
117 libmesh_error_msg(
"Error! Zero area_swept_function specified!" <<
120 this->set_parameter(aoa_function, input,
124 if (aoa_function.expression() ==
"00000")
125 libmesh_error_msg(
"Error! No angle-of-attack specified!" <<
128 this->set_parameter(torque_function, input,
132 if (torque_function.expression() ==
"0")
133 std::cout <<
"Warning! Zero torque function specified!" << std::endl;
136 (this->moment_of_inertia, input,
140 if (!moment_of_inertia)
142 "Error! Zero AveragedTurbine moment of inertia specified!" <<
146 (this->initial_speed, input,
150 this->_fan_speed_var_name = input(
"Physics/VariableNames/fan_speed",
157 (
const libMesh::Point& point,
158 const libMesh::Real time,
159 const libMesh::NumberVectorValue& U,
161 libMesh::NumberVectorValue& U_B_1,
162 libMesh::NumberVectorValue& F,
163 libMesh::NumberTensorValue *dFdU,
164 libMesh::NumberVectorValue *dFds)
167 libMesh::DenseVector<libMesh::Number> output_vec(3);
169 base_velocity_function(point, time, output_vec);
171 U_B_1(0) = output_vec(0);
172 U_B_1(1) = output_vec(1);
173 U_B_1(2) = output_vec(2);
175 const libMesh::NumberVectorValue U_B = U_B_1 * s;
177 const libMesh::Number U_B_size = U_B.size();
184 const libMesh::NumberVectorValue N_B =
185 libMesh::NumberVectorValue(U_B/U_B_size);
187 local_vertical_function(point, time, output_vec);
190 const libMesh::NumberVectorValue N_V(output_vec(0),
196 const libMesh::NumberVectorValue N_R = N_B.cross(N_V);
199 const libMesh::NumberVectorValue U_P = U - (U*N_R)*N_R - U_B;
201 const libMesh::Number U_P_size = U_P.size();
210 const libMesh::NumberVectorValue N_drag =
211 libMesh::NumberVectorValue(-U_P/U_P_size);
214 const libMesh::NumberVectorValue N_lift = N_drag.cross(N_R);
217 const libMesh::Number u_fwd = -(U_P * N_B);
220 const libMesh::Number u_up = U_P * N_V;
224 libmesh_assert (u_up || u_fwd);
227 const libMesh::Number part_angle = std::atan2(u_up, u_fwd);
230 const libMesh::Number angle = part_angle +
231 aoa_function(point, time);
233 const libMesh::Number C_lift = lift_function(point, angle);
234 const libMesh::Number C_drag = drag_function(point, angle);
236 const libMesh::Number chord = chord_function(point, time);
237 const libMesh::Number area = area_swept_function(point, time);
239 const libMesh::Number v_sq = U_P*U_P;
241 const libMesh::Number LDfactor = 0.5 * this->_rho * v_sq * chord / area;
242 const libMesh::Number lift = C_lift * LDfactor;
243 const libMesh::Number drag = C_drag * LDfactor;
246 F = lift * N_lift + drag * N_drag;
250 const libMesh::NumberVectorValue LDderivfactor =
251 (N_lift*C_lift+N_drag*C_drag) *
252 this->_rho * chord / area;
254 const libMesh::Number sfactor = -(U_P*U_B_1);
256 (*dFds) = LDderivfactor * sfactor;
258 for (
unsigned int i=0; i != 3; ++i)
259 for (
unsigned int j=0; j != 3; ++j)
260 (*dFdU)(i,j) = LDderivfactor(i) * U_P(j);
virtual void set_time_evolving_vars(libMesh::FEMSystem *system)
Sets turbine_speed and velocity variables to be time-evolving.
virtual void init_variables(libMesh::FEMSystem *system)
Initialization of variables.
Physics class for Incompressible Navier-Stokes.
const std::string fan_speed_var_name_default
fan speed
INSTANTIATE_INC_NS_SUBCLASS(AveragedTurbineBase)
virtual void init_variables(libMesh::FEMSystem *system)
Initialization of Navier-Stokes variables.
static PhysicsName averaged_turbine()
virtual void set_time_evolving_vars(libMesh::FEMSystem *system)
Sets velocity variables to be time-evolving.
bool compute_force(const libMesh::Point &point, const libMesh::Real time, const libMesh::NumberVectorValue &U, libMesh::Number s, libMesh::NumberVectorValue &U_B_1, libMesh::NumberVectorValue &F, libMesh::NumberTensorValue *dFdU=NULL, libMesh::NumberVectorValue *dFds=NULL)
void read_input_options(const GetPot &input)
Read options from GetPot input file.