adjointShapeOptimizationFoam.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2017 OpenFOAM Foundation
9-------------------------------------------------------------------------------
10License
11 This file is part of OpenFOAM.
12
13 OpenFOAM is free software: you can redistribute it and/or modify it
14 under the terms of the GNU General Public License as published by
15 the Free Software Foundation, either version 3 of the License, or
16 (at your option) any later version.
17
18 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 for more details.
22
23 You should have received a copy of the GNU General Public License
24 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25
26Application
27 ajointShapeOptimizationFoam
28
29Group
30 grpIncompressibleSolvers
31
32Description
33 Steady-state solver for incompressible, turbulent flow of non-Newtonian
34 fluids with optimisation of duct shape by applying "blockage" in regions
35 causing pressure loss as estimated using an adjoint formulation.
36
37 References:
38 \verbatim
39 "Implementation of a continuous adjoint for topology optimization of
40 ducted flows"
41 C. Othmer,
42 E. de Villiers,
43 H.G. Weller
44 AIAA-2007-3947
45 http://pdf.aiaa.org/preview/CDReadyMCFD07_1379/PV2007_3947.pdf
46 \endverbatim
47
48 Note that this solver optimises for total pressure loss whereas the
49 above paper describes the method for optimising power-loss.
50
51\*---------------------------------------------------------------------------*/
52
53#include "fvCFD.H"
56#include "simpleControl.H"
57#include "fvOptions.H"
58
59template<class Type>
60void zeroCells
61(
62 GeometricField<Type, fvPatchField, volMesh>& vf,
63 const labelList& cells
64)
65{
66 forAll(cells, i)
67 {
68 vf[cells[i]] = Zero;
69 }
70}
71
72
73// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
74
75int main(int argc, char *argv[])
76{
77 argList::addNote
78 (
79 "Steady-state solver for incompressible, turbulent flow"
80 " of non-Newtonian fluids with duct shape optimisation"
81 " by applying 'blockage' in regions causing pressure loss"
82 );
83
84 #include "postProcess.H"
85
86 #include "addCheckCaseOptions.H"
87 #include "setRootCaseLists.H"
88 #include "createTime.H"
89 #include "createMesh.H"
90 #include "createControl.H"
91 #include "createFields.H"
92 #include "initContinuityErrs.H"
94
95 turbulence->validate();
96
97 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
98
99 Info<< "\nStarting time loop\n" << endl;
100
101 while (simple.loop())
102 {
103 Info<< "Time = " << runTime.timeName() << nl << endl;
104
105 //alpha +=
106 // mesh.relaxationFactor("alpha")
107 // *(lambda*max(Ua & U, zeroSensitivity) - alpha);
108 alpha +=
109 mesh.fieldRelaxationFactor("alpha")
110 *(min(max(alpha + lambda*(Ua & U), zeroAlpha), alphaMax) - alpha);
111
113 //zeroCells(alpha, outletCells);
114
115 // Pressure-velocity SIMPLE corrector
116 {
117 // Momentum predictor
118
119 tmp<fvVectorMatrix> tUEqn
120 (
121 fvm::div(phi, U)
122 + turbulence->divDevReff(U)
123 + fvm::Sp(alpha, U)
124 ==
125 fvOptions(U)
126 );
127 fvVectorMatrix& UEqn = tUEqn.ref();
128
129 UEqn.relax();
130
131 fvOptions.constrain(UEqn);
132
133 solve(UEqn == -fvc::grad(p));
134
135 fvOptions.correct(U);
136
137 volScalarField rAU(1.0/UEqn.A());
139 tUEqn.clear();
140 surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
141 adjustPhi(phiHbyA, U, p);
142
143 // Update the pressure BCs to ensure flux consistency
145
146 // Non-orthogonal pressure corrector loop
147 while (simple.correctNonOrthogonal())
148 {
149 fvScalarMatrix pEqn
150 (
151 fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
152 );
153
154 pEqn.setReference(pRefCell, pRefValue);
155 pEqn.solve();
156
157 if (simple.finalNonOrthogonalIter())
158 {
159 phi = phiHbyA - pEqn.flux();
160 }
161 }
162
163 #include "continuityErrs.H"
164
165 // Explicitly relax pressure for momentum corrector
166 p.relax();
167
168 // Momentum corrector
169 U = HbyA - rAU*fvc::grad(p);
170 U.correctBoundaryConditions();
171 fvOptions.correct(U);
172 }
173
174 // Adjoint Pressure-velocity SIMPLE corrector
175 {
176 // Adjoint Momentum predictor
177
178 volVectorField adjointTransposeConvection((fvc::grad(Ua) & U));
179 //volVectorField adjointTransposeConvection
180 //(
181 // fvc::reconstruct
182 // (
183 // mesh.magSf()*fvc::dotInterpolate(fvc::snGrad(Ua), U)
184 // )
185 //);
186
187 zeroCells(adjointTransposeConvection, inletCells);
188
189 tmp<fvVectorMatrix> tUaEqn
190 (
191 fvm::div(-phi, Ua)
192 - adjointTransposeConvection
193 + turbulence->divDevReff(Ua)
194 + fvm::Sp(alpha, Ua)
195 ==
196 fvOptions(Ua)
197 );
198 fvVectorMatrix& UaEqn = tUaEqn.ref();
199
200 UaEqn.relax();
201
202 fvOptions.constrain(UaEqn);
203
204 solve(UaEqn == -fvc::grad(pa));
205
206 fvOptions.correct(Ua);
207
208 volScalarField rAUa(1.0/UaEqn.A());
209 volVectorField HbyAa("HbyAa", Ua);
210 HbyAa = rAUa*UaEqn.H();
211 tUaEqn.clear();
212 surfaceScalarField phiHbyAa("phiHbyAa", fvc::flux(HbyAa));
213 adjustPhi(phiHbyAa, Ua, pa);
214
215 // Non-orthogonal pressure corrector loop
216 while (simple.correctNonOrthogonal())
217 {
219 (
220 fvm::laplacian(rAUa, pa) == fvc::div(phiHbyAa)
221 );
222
223 paEqn.setReference(paRefCell, paRefValue);
224 paEqn.solve();
225
226 if (simple.finalNonOrthogonalIter())
227 {
228 phia = phiHbyAa - paEqn.flux();
229 }
230 }
231
232 #include "adjointContinuityErrs.H"
233
234 // Explicitly relax pressure for adjoint momentum corrector
235 pa.relax();
236
237 // Adjoint momentum corrector
238 Ua = HbyAa - rAUa*fvc::grad(pa);
239 Ua.correctBoundaryConditions();
240 fvOptions.correct(Ua);
241 }
242
243 laminarTransport.correct();
244 turbulence->correct();
245
246 runTime.write();
247
248 runTime.printExecutionTime(Info);
249 }
250
251 Info<< "End\n" << endl;
252
253 return 0;
254}
255
256
257// ************************************************************************* //
Y[inertIndex] max(0.0)
Required Classes.
Calculates and prints the continuity errors.
fv::options & fvOptions
surfaceScalarField & phi
const scalar pRefValue
const label pRefCell
U
Definition: pEqn.H:72
volScalarField & p
tmp< fvVectorMatrix > tUEqn(fvm::ddt(rho, U)+fvm::div(phi, U)+MRF.DDt(rho, U)+turbulence->divDevRhoReff(U)==fvOptions(rho, U))
fvVectorMatrix & UEqn
Definition: UEqn.H:13
phiHbyA
Definition: pcEqn.H:73
HbyA
Definition: pcEqn.H:74
dynamicFvMesh & mesh
engineTime & runTime
Required Variables.
compressible::turbulenceModel & turbulence
const cellShapeList & cells
adjustPhi(phiHbyA, U, p_rgh)
Declare and initialise the cumulative ddjoint continuity error.
tmp< volScalarField > rAU
Definition: initCorrectPhi.H:1
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:83
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:45
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:82
messageStream Info
Information stream (stdout output on master, null elsewhere)
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
fvMatrix< vector > fvVectorMatrix
Definition: fvMatricesFwd.H:46
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
fvScalarMatrix paEqn(fvm::d2dt2(pa) - sqr(c0) *fvc::laplacian(pa))
Execute application functionObjects to post-process existing results.
const dictionary & simple
volScalarField & alpha
CEqn solve()
singlePhaseTransportModel laminarTransport(U, phi)
dimensionedScalar zeroAlpha(dimless/dimTime, Zero)
const labelList & inletCells
Definition: createFields.H:106
dimensionedScalar lambda("lambda", dimTime/sqr(dimLength), laminarTransport)
zeroCells(alpha, inletCells)
dimensionedScalar alphaMax("alphaMax", dimless/dimTime, laminarTransport)
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333