From 7f75ace2ff6e82dfa699b9e8775cf8947f54d57f Mon Sep 17 00:00:00 2001 From: Josh Anibal Date: Sat, 15 Jun 2024 08:36:24 -0700 Subject: [PATCH] Debug update of code to 3.40 (#35) * added new src 1 * added new src 2 * added new src 2 * added new src 2 * added new src * added atime.f but it is not used' * added unitset' * added new src amode * added new src avl.f * added more variables but reduced size * updated amake * added remaining changes to AVL.inc * fixed *DIR error in exec newton system * bump version * updated ad files --------- Co-authored-by: josh --- meson.build | 2 +- pyavl/pyAVL.py | 15 +- pyproject.toml | 2 +- src/ad_src/forward_ad_src/aero_d.f | 66 ++- src/ad_src/forward_ad_src/aic_d.f | 2 + src/ad_src/forward_ad_src/amake_d.f | 126 ++++- src/ad_src/forward_ad_src/atpforc_d.f | 71 +-- src/ad_src/reverse_ad_src/aero_b.f | 73 +-- src/ad_src/reverse_ad_src/aic_b.f | 2 + src/ad_src/reverse_ad_src/amake_b.f | 52 +- src/ad_src/reverse_ad_src/atpforc_b.f | 55 ++- src/aero.f | 53 ++- src/aic.f | 5 +- src/ainput.f | 410 +++++++++++----- src/amake.f | 111 ++++- src/amode.f | 93 ++-- src/aoml.f | 470 ++++++++++++++++++ src/aoper.f | 134 ++++-- src/aoutmrf.f | 587 +++++++++++++++++++++++ src/aoutput.f | 573 ++++++++++++++++++---- src/atime.f | 656 ++++++++++++++++++++++++++ src/atpforc.f | 43 +- src/atrim.f | 134 +++--- src/avl.f | 129 ++--- src/build/fileList | 2 + src/f2py/libavl.pyf | 8 +- src/includes/AVL.INC | 165 ++++--- src/includes/AVLPLT.INC | 1 + src/unitset.f | 486 +++++++++++++++++++ tests/test_analysis.py | 8 +- tests/test_io.py | 5 +- 31 files changed, 3871 insertions(+), 668 deletions(-) create mode 100644 src/aoml.f create mode 100644 src/aoutmrf.f create mode 100644 src/atime.f create mode 100644 src/unitset.f diff --git a/meson.build b/meson.build index a8cc997..61e4493 100644 --- a/meson.build +++ b/meson.build @@ -1,7 +1,7 @@ project( 'pyavl', 'c', - version: '1.6.4', + version: '1.7.0', license: 'GPL-2.0', meson_version: '>= 0.64.0', default_options: [ diff --git a/pyavl/pyAVL.py b/pyavl/pyAVL.py index 607464f..16b62d1 100644 --- a/pyavl/pyAVL.py +++ b/pyavl/pyAVL.py @@ -754,8 +754,9 @@ def executeRun(self): warnings.warn("executeRun is deprecated, use execute_run instead") self.execute_run() - def execute_run(self): + def execute_run(self, tol=0.00002): # run the analysis (equivalent to the avl command `x` in the oper menu) + self.set_avl_fort_arr('CASE_R', 'EXEC_TOL', tol) self.avl.oper() def CLSweep(self, start_CL, end_CL, increment=0.1): @@ -1151,8 +1152,16 @@ def __write_data(key_list, newline: bool = True): fid.write(f" {data['claf'][idx_sec]}\n") if (data["clcdsec"][idx_sec] != 0.0).any(): - fid.write(" CLCD\n") - fid.write(f" {data['clcd'][idx_sec]}\n") + fid.write(" CDCL\n") + fid.write( + f" {data['clcdsec'][idx_sec, 0]:.6f} " + f" {data['clcdsec'][idx_sec, 1]:.6f} " + f" {data['clcdsec'][idx_sec, 2]:.6f} " + f" {data['clcdsec'][idx_sec, 3]:.6f} " + f" {data['clcdsec'][idx_sec, 4]:.6f} " + f" {data['clcdsec'][idx_sec, 5]:.6f}\n" + ) + # check for control surfaces diff --git a/pyproject.toml b/pyproject.toml index 6ec6cca..3733a7a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -19,7 +19,7 @@ dependencies = [ "numpy>=1.19", ] readme = "README.md" -version = "1.6.4" # this automatically updates __init__.py +version = "1.7.0" # this automatically updates __init__.py [tool.cibuildwheel] diff --git a/src/ad_src/forward_ad_src/aero_d.f b/src/ad_src/forward_ad_src/aero_d.f index d7c5078..f684f30 100644 --- a/src/ad_src/forward_ad_src/aero_d.f +++ b/src/ad_src/forward_ad_src/aero_d.f @@ -1343,9 +1343,15 @@ SUBROUTINE SFFORC_D() C C******************************************************************* C -C...Store strip X,Y,Z body axes forces -C (these are normalized by strip area and moments are referred to -C c/4 point and are normalized by strip chord and area) +C--- At this point the forces are accumulated for the strip in body axes, +C referenced to strip 1/4 chord and normalized by strip area and chord +C CFX, CFY, CFZ ! body axes forces +C CMX, CMY, CMZ ! body axes moments about c/4 point +C CNC ! strip spanloading CN*chord +C CDV_LSTRP ! strip viscous drag in stability axes +C +C...Store strip X,Y,Z body axes forces and moments +C referred to c/4 and strip area and chord cf_strp(1, j) = cfx cf_strp(2, j) = cfy cf_strp(3, j) = cfz @@ -1353,19 +1359,21 @@ SUBROUTINE SFFORC_D() cm_strp(2, j) = cmy cm_strp(3, j) = cmz C -C...Transform strip body axes forces into stability axes - cdstrp_diff(j) = cosa*cfx_diff + cfx*cosa_diff + sina*cfz_diff + - + cfz*sina_diff - cdstrp(j) = cfx*cosa + cfz*sina - clstrp_diff(j) = cosa*cfz_diff + cfz*cosa_diff - sina*cfx_diff - - + cfx*sina_diff - clstrp(j) = -(cfx*sina) + cfz*cosa +C...Strip body axes forces, referred to strip area and chord cxstrp_diff(j) = cfx_diff cxstrp(j) = cfx cystrp_diff(j) = cfy_diff cystrp(j) = cfy czstrp_diff(j) = cfz_diff czstrp(j) = cfz +C...Transform strip body axes forces into stability axes, +C referred to strip area and chord + cdstrp_diff(j) = cosa*cfx_diff + cfx*cosa_diff + sina*cfz_diff + + + cfz*sina_diff + cdstrp(j) = cfx*cosa + cfz*sina + clstrp_diff(j) = cosa*cfz_diff + cfz*cosa_diff - sina*cfx_diff - + + cfx*sina_diff + clstrp(j) = -(cfx*sina) + cfz*cosa C cdst_a(j) = -(cfx*sina) + cfz*cosa clst_a(j) = -(cfx*cosa) - cfz*sina @@ -1401,14 +1409,15 @@ SUBROUTINE SFFORC_D() czst_g(j, n) = cfz_g(n) ENDDO C -C... Set strip moments about the overall moment reference point XYZREF -C (still normalized by strip area and chord) +C------ vector from chord c/4 reference point to case reference point XYZREF r_diff(1) = xr_diff - xyzref_diff(1) r(1) = xr - xyzref(1) r_diff(2) = yr_diff - xyzref_diff(2) r(2) = yr - xyzref(2) r_diff(3) = zr_diff - xyzref_diff(3) r(3) = zr - xyzref(3) +C... Strip moments in body axes about the case moment reference point XYZREF +C normalized by strip area and chord temp = (cfz*r(2)-cfy*r(3))/cr crstrp_diff(j) = cmx_diff + (r(2)*cfz_diff+cfz*r_diff(2)-r(3)* + cfy_diff-cfy*r_diff(3)-temp*cr_diff)/cr @@ -1452,9 +1461,8 @@ SUBROUTINE SFFORC_D() cnst_g(j, n) = cmz_g(n) + (cfy_g(n)*r(1)-cfx_g(n)*r(2))/cr ENDDO C -C...Take components of X,Y,Z forces in local strip axes -C (axial/normal and lift/drag) -C in plane normal to (possibly dihedralled) strip +C...Components of X,Y,Z forces in local strip axes +C axial/normal forces and lift/drag in plane normal to dihedral of strip cl_lstrp(j) = ulift(1)*cfx + ulift(2)*cfy + ulift(3)*cfz cd_lstrp(j) = udrag(1)*cfx + udrag(2)*cfy + udrag(3)*cfz caxlstrp(j) = cfx @@ -1468,6 +1476,7 @@ SUBROUTINE SFFORC_D() rrot(2) = ysref(j) - xyzref(2) rrot_diff(3) = zsref_diff(j) - xyzref_diff(3) rrot(3) = zsref(j) - xyzref(3) +C print *,"WROT ",WROT C C------ set total effective velocity = freestream + rotation DO ii1=1,3 @@ -1718,6 +1727,7 @@ SUBROUTINE SFFORC_D() enave(2) = enave(2) + sr*ensy(j) enave(3) = enave(3) + sr*ensz(j) C +C--- Surface lift and drag referenced to case SREF, CREF, BREF temp0 = sr/sref cdsurf_diff(is) = cdsurf_diff(is) + temp0*cdstrp_diff(j) + + cdstrp(j)*(sr_diff-temp0*sref_diff)/sref @@ -1726,7 +1736,7 @@ SUBROUTINE SFFORC_D() clsurf_diff(is) = clsurf_diff(is) + temp0*clstrp_diff(j) + + clstrp(j)*(sr_diff-temp0*sref_diff)/sref clsurf(is) = clsurf(is) + clstrp(j)*temp0 -C +C--- Surface body axes forces referenced to case SREF, CREF, BREF temp0 = sr/sref cxsurf_diff(is) = cxsurf_diff(is) + temp0*cxstrp_diff(j) + + cxstrp(j)*(sr_diff-temp0*sref_diff)/sref @@ -1739,7 +1749,7 @@ SUBROUTINE SFFORC_D() czsurf_diff(is) = czsurf_diff(is) + temp0*czstrp_diff(j) + + czstrp(j)*(sr_diff-temp0*sref_diff)/sref czsurf(is) = czsurf(is) + czstrp(j)*temp0 -C +C--- Surface body axes moments referenced to case SREF, CREF, BREF about XYZREF temp0 = crstrp(j)*sr*cr/(sref*bref) crsurf_diff(is) = crsurf_diff(is) + (sr*cr*crstrp_diff(j)+ + crstrp(j)*(cr*sr_diff+sr*cr_diff)-temp0*(bref*sref_diff+sref @@ -1757,6 +1767,7 @@ SUBROUTINE SFFORC_D() cnsurf(is) = cnsurf(is) + temp0 C C--- Bug fix, HHY/S.Allmaras +C--- Surface viscous drag referenced to case SREF, CREF, BREF temp0 = sr/sref cdvsurf_diff(is) = cdvsurf_diff(is) + temp0*cdv_lstrp_diff(j) + + cdv_lstrp(j)*(sr_diff-temp0*sref_diff)/sref @@ -1899,32 +1910,35 @@ SUBROUTINE SFFORC_D() END IF cl_srf(is) = DOT(ulift, cf_srf(1, is)) cd_srf(is) = DOT(udrag, cf_srf(1, is)) -C--- Surface hinge moments defined by surface LE moment about hinge vector +C +C--- Surface hinge moments defined by surface LE moment about hinge vector Ccc CMLE_SRF(IS) = DOT(CM_SRF(1,IS),VHINGE(1,IS)) C C C------------------------------------------------- IF (lfload(is)) THEN -C------- Total forces summed from surface forces... -C- normalized to configuration reference quantities - cdtot_diff = cdtot_diff + cdsurf_diff(is) - cdtot = cdtot + cdsurf(is) - cltot_diff = cltot_diff + clsurf_diff(is) - cltot = cltot + clsurf(is) +C--- Total forces summed from surface forces +C normalized to case reference quantities SREF, CREF, BREF cxtot_diff = cxtot_diff + cxsurf_diff(is) cxtot = cxtot + cxsurf(is) cytot_diff = cytot_diff + cysurf_diff(is) cytot = cytot + cysurf(is) cztot_diff = cztot_diff + czsurf_diff(is) cztot = cztot + czsurf(is) + cdtot_diff = cdtot_diff + cdsurf_diff(is) + cdtot = cdtot + cdsurf(is) + cltot_diff = cltot_diff + clsurf_diff(is) + cltot = cltot + clsurf(is) + cdvtot_diff = cdvtot_diff + cdvsurf_diff(is) + cdvtot = cdvtot + cdvsurf(is) +C--- Total body axes moments about XYZREF summed from surface moments +C normalized to case reference quantities SREF, CREF, BREF crtot_diff = crtot_diff + crsurf_diff(is) crtot = crtot + crsurf(is) cmtot_diff = cmtot_diff + cmsurf_diff(is) cmtot = cmtot + cmsurf(is) cntot_diff = cntot_diff + cnsurf_diff(is) cntot = cntot + cnsurf(is) - cdvtot_diff = cdvtot_diff + cdvsurf_diff(is) - cdvtot = cdvtot + cdvsurf(is) C cdtot_a = cdtot_a + cds_a(is) cltot_a = cltot_a + cls_a(is) diff --git a/src/ad_src/forward_ad_src/aic_d.f b/src/ad_src/forward_ad_src/aic_d.f index c713b37..fdc56a2 100644 --- a/src/ad_src/forward_ad_src/aic_d.f +++ b/src/ad_src/forward_ad_src/aic_d.f @@ -370,6 +370,8 @@ SUBROUTINE VORVELC_D(x, x_diff, y, y_diff, z, z_diff, lbound, x1, + rcore, rcore_diff) C---------------------------------------------------------- C Same as VORVEL, with finite core radius +C Uses Scully (also Burnham-Hallock) core model +C Vtan = Gam/2*pi . r/(r^2 +rcore^2) C---------------------------------------------------------- LOGICAL lbound C diff --git a/src/ad_src/forward_ad_src/amake_d.f b/src/ad_src/forward_ad_src/amake_d.f index 8346757..5a683b4 100644 --- a/src/ad_src/forward_ad_src/amake_d.f +++ b/src/ad_src/forward_ad_src/amake_d.f @@ -214,6 +214,7 @@ SUBROUTINE MAKESURF_D(isurf) REAL ypt1_diff REAL yscale REAL yscale_diff + INTEGER ii REAL width REAL width_diff REAL chordl @@ -293,6 +294,8 @@ SUBROUTINE MAKESURF_D(isurf) REAL fracte_diff INTRINSIC MAX INTRINSIC MIN + REAL zl + REAL zu REAL sum REAL wtot INTEGER jj @@ -360,7 +363,7 @@ SUBROUTINE MAKESURF_D(isurf) idx_strip = jfrst(isurf) C C----------------------------------------------------------------- -C---- section arc lengths of wing trace in y-z plane +C---- Arc length positions of sections in wing trace in y-z plane yzlen(1) = 0. DO ii1=1,ksmax yzlen_diff(ii1) = 0.D0 @@ -446,7 +449,13 @@ SUBROUTINE MAKESURF_D(isurf) END IF ELSE C -C----- set spanwise spacing using overall parameters NVS(ISURF), SSPACE +C +C----- Otherwise, set spanwise spacing using the SURFACE spanwise +C parameters NVS, SSPACE +C +C This spanwise spacing is modified (fudged) to align vortex edges +C with SECTIONs as defined. This allows CONTROLs to be defined +C without bridging vortex strips C nspace = 2*nvs(isurf) + 1 IF (nspace .GT. kpmax) THEN @@ -497,7 +506,7 @@ SUBROUTINE MAKESURF_D(isurf) iptloc(1) = 1 iptloc(nsec(isurf)) = npt C -C----- fudge Glauert angles to make nodes match up exactly with interior sections +C----- fudge spacing array to make nodes match up exactly with interior sections DO isec=2,nsec(isurf)-1 ipt1 = iptloc(isec-1) ipt2 = iptloc(isec) @@ -505,6 +514,7 @@ SUBROUTINE MAKESURF_D(isurf) GOTO 110 ELSE C +C----- fudge spacing to this section so that nodes match up exactly with section ypt1_diff = ypt_diff(ipt1) ypt1 = ypt(ipt1) temp = (yzlen(isec)-yzlen(isec-1))/(ypt(ipt2)-ypt(ipt1)) @@ -522,12 +532,14 @@ SUBROUTINE MAKESURF_D(isurf) ycp(ivs) = yzlen(isec-1) + yscale*(ycp(ivs)-ypt1) ENDDO C +C----- check for unique spacing node for next section, if not we need more nodes ipt1 = iptloc(isec) ipt2 = iptloc(isec+1) IF (ipt1 .EQ. ipt2) THEN GOTO 120 ELSE C +C----- fudge spacing to this section so that nodes match up exactly with section ypt1_diff = ypt_diff(ipt1) ypt1 = ypt(ipt1) temp = (ypt(ipt2)-yzlen(isec))/(ypt(ipt2)-ypt(ipt1)) @@ -558,12 +570,25 @@ SUBROUTINE MAKESURF_D(isurf) END IF END IF C +Cc#ifdef USE_CPOML +C... store section counters + 130 IF (isurf .EQ. 1) THEN + icntfrst(isurf) = 1 + ELSE + icntfrst(isurf) = icntfrst(isurf-1) + ncntsec(isurf-1) + END IF + ncntsec(isurf) = nsec(isurf) + DO isec=1,nsec(isurf) + ii = icntfrst(isurf) + (isec-1) + icntsec(ii) = iptloc(isec) + ENDDO +Cc#endif C C C==================================================== C---- define strips between input sections C - 130 nj(isurf) = 0 + nj(isurf) = 0 C IF (ncontrol .GT. ndmax) THEN WRITE(*, *) @@ -823,6 +848,15 @@ SUBROUTINE MAKESURF_D(isurf) tante(idx_strip) = (xyzler(1)+chordr-xyzlel(1)-chordl)/ + width C +Cc#ifdef USE_CPOML + chsin = chsinl + f1*(chsinr-chsinl) + chcos = chcosl + f1*(chcosr-chcosl) + ainc1(idx_strip) = ATAN2(chsin, chcos) + chsin = chsinl + f2*(chsinr-chsinl) + chcos = chcosl + f2*(chcosr-chcosl) + ainc2(idx_strip) = ATAN2(chsin, chcos) +C +Cc#endif chsin_diff = chsinl_diff + (chsinr-chsinl)*fc_diff + fc*( + chsinr_diff-chsinl_diff) chsin = chsinl + fc*(chsinr-chsinl) @@ -1032,6 +1066,8 @@ SUBROUTINE MAKESURF_D(isurf) C-------- go over vortices in this strip idx_vor = ijfrst(idx_strip) C NVOR = NVOR + 1 +C change all NVOR indices into idx_vor +C change all NSTRIP indices into idx_strip DO ivc=1,nvc(isurf) C rv1_diff(1, idx_vor) = rle1_diff(1, idx_strip) + xvr(ivc @@ -1172,9 +1208,40 @@ SUBROUTINE MAKESURF_D(isurf) C C---------- TE control point used only if surface sheds a wake lvnc(idx_vor) = lfwake(isurf) +C +Cc#ifdef USE_CPOML +C... nodal grid associated with vortex strip (aft-panel nodes) +C... NOTE: airfoil in plane of wing, but not rotated perpendicular to dihedral; +C... retained in (x,z) plane at this point + CALL AKIMA(xlasec(1, isec, isurf), zlasec(1, isec, isurf + + ), nsl, xpt(ivc+1), zl, dsdx) + CALL AKIMA(xuasec(1, isec, isurf), zuasec(1, isec, isurf + + ), nsl, xpt(ivc+1), zu, dsdx) +C + xyn1(1, idx_vor) = rle1(1, idx_strip) + xpt(ivc+1)* + + chord1(idx_strip) + xyn1(2, idx_vor) = rle1(2, idx_strip) + zlon1(idx_vor) = rle1(3, idx_strip) + zl*chord1( + + idx_strip) + zupn1(idx_vor) = rle1(3, idx_strip) + zu*chord1( + + idx_strip) +C + CALL AKIMA(xlasec(1, isec+1, isurf), zlasec(1, isec+1, + + isurf), nsl, xpt(ivc+1), zl, dsdx) + CALL AKIMA(xuasec(1, isec+1, isurf), zuasec(1, isec+1, + + isurf), nsl, xpt(ivc+1), zu, dsdx) +C + xyn2(1, idx_vor) = rle2(1, idx_strip) + xpt(ivc+1)* + + chord2(idx_strip) + xyn2(2, idx_vor) = rle2(2, idx_strip) + zlon2(idx_vor) = rle2(3, idx_strip) + zl*chord2( + + idx_strip) + zupn2(idx_vor) = rle2(3, idx_strip) + zu*chord2( + + idx_strip) +C +Cc#endif idx_vor = idx_vor + 1 ENDDO -C C idx_strip = idx_strip + 1 ENDDO @@ -1233,6 +1300,9 @@ SUBROUTINE SDUPL_D(nn, ypt, msg) INTEGER klen INTRINSIC LEN INTEGER k + INTEGER isec + INTEGER idup + INTEGER iorg REAL yoff INTEGER idx_strip INTEGER ivs @@ -1250,7 +1320,8 @@ SUBROUTINE SDUPL_D(nn, ypt, msg) C nni = nn + 1 IF (nni .GT. nfmax) THEN - WRITE(*, *) 'SDUPL: Surface array overflow. Increase NFMAX.' + WRITE(*, *) 'SDUPL: Surface array overflow. Increase NFMAX', + + ' currently ', nfmax STOP ELSE C @@ -1272,6 +1343,7 @@ SUBROUTINE SDUPL_D(nn, ypt, msg) lfwake(nni) = lfwake(nn) lfalbe(nni) = lfalbe(nn) lfload(nni) = lfload(nn) + lrange(nni) = lrange(nn) C IFRST(NNI) = NVOR + 1 lsurfspacing(nni) = lsurfspacing(nn) C @@ -1293,6 +1365,16 @@ SUBROUTINE SDUPL_D(nn, ypt, msg) C C--- Image flag reversed (set to -IMAGS) for imaged surfaces imags(nni) = -imags(nn) +C +Cc#ifdef USE_CPOML + icntfrst(nni) = icntfrst(nn) + ncntsec(nn) + ncntsec(nni) = ncntsec(nn) + DO isec=1,ncntsec(nni) + idup = icntfrst(nni) + (isec-1) + iorg = icntfrst(nn) + (isec-1) + icntsec(idup) = icntsec(iorg) + ENDDO +Cc#endif C yoff = 2.0*ypt C @@ -1303,7 +1385,7 @@ SUBROUTINE SDUPL_D(nn, ypt, msg) C NSTRIP = NSTRIP + 1 DO ivs=1,nvs(nni) IF (idx_strip .GT. nsmax) THEN - GOTO 110 + GOTO 100 ELSE C jji = jfrst(nni) + ivs - 1 @@ -1338,6 +1420,11 @@ SUBROUTINE SDUPL_D(nn, ypt, msg) ainc_diff(jji) = ainc_diff(jj) ainc(jji) = ainc(jj) C +Cc#ifdef USE_CPOML + ainc1(jji) = ainc2(jj) + ainc2(jji) = ainc1(jj) +C +Cc#endif nsurfs(idx_strip) = nni C DO n=1,ndesign @@ -1376,7 +1463,7 @@ SUBROUTINE SDUPL_D(nn, ypt, msg) C NVOR = NVOR + 1 DO ivc=1,nvc(nni) IF (idx_vor .GT. nvmax) THEN - GOTO 120 + GOTO 110 ELSE C iii = ijfrst(jji) + ivc - 1 @@ -1423,23 +1510,38 @@ SUBROUTINE SDUPL_D(nn, ypt, msg) dcontrol_diff(iii, n) = -(rsgn*dcontrol_diff(ii, n)) dcontrol(iii, n) = -(dcontrol(ii, n)*rsgn) ENDDO +C +Cc#ifdef USE_CPOML +C... nodal grid associated with vortex strip + xyn1(1, iii) = xyn2(1, ii) + xyn1(2, iii) = -xyn2(2, ii) + yoff + xyn2(1, iii) = xyn1(1, ii) + xyn2(2, iii) = -xyn1(2, ii) + yoff +C + zlon1(iii) = zlon2(ii) + zupn1(iii) = zupn2(ii) + zlon2(iii) = zlon1(ii) + zupn2(iii) = zupn1(ii) +Cc#endif idx_vor = idx_vor + 1 END IF ENDDO +C + idx_strip = idx_strip + 1 END IF ENDDO C C -C -C C nstrip = nstrip + nj(nni) nvor = nvor + nk(nni)*nj(nni) C RETURN - 110 WRITE(*, *) 'SDUPL: Strip array overflow. Increase NSMAX.' + 100 WRITE(*, *) 'SDUPL: Strip array overflow. Increase NSMAX', + + ' currently ', nsmax STOP - 120 WRITE(*, *) 'SDUPL: Vortex array overflow. Increase NVMAX.' + 110 WRITE(*, *) 'SDUPL: Vortex array overflow. Increase NVMAX', + + ' currently ', nvmax STOP END IF END diff --git a/src/ad_src/forward_ad_src/atpforc_d.f b/src/ad_src/forward_ad_src/atpforc_d.f index a9a9743..a1885d1 100644 --- a/src/ad_src/forward_ad_src/atpforc_d.f +++ b/src/ad_src/forward_ad_src/atpforc_d.f @@ -79,6 +79,7 @@ SUBROUTINE TPFORC_D() REAL rsq1_diff REAL rsq2 REAL rsq2_diff + INTEGER isurf REAL ar REAL ar_diff REAL spanef_cl @@ -417,40 +418,46 @@ SUBROUTINE TPFORC_D() ENDDO C C -C...Trefftz-plane drag is kinetic energy in crossflow dwwake(jc) = -(ny*vy+nz*vz) C - temp1 = dyt/sref - clff_diff = clff_diff + 2.0*(temp1*gams_diff(jc)+gams(jc)*( - + dyt_diff-temp1*sref_diff)/sref) - clff = clff + 2.0*(gams(jc)*temp1) - temp1 = dzt/sref - cyff_diff = cyff_diff - 2.0*(temp1*gams_diff(jc)+gams(jc)*( - + dzt_diff-temp1*sref_diff)/sref) - cyff = cyff - 2.0*(gams(jc)*temp1) - temp0 = dzt*vy - dyt*vz - temp1 = gams(jc)/sref - cdff_diff = cdff_diff + temp0*(gams_diff(jc)-temp1*sref_diff)/ - + sref + temp1*(vy*dzt_diff+dzt*vy_diff-vz*dyt_diff-dyt*vz_diff) - cdff = cdff + temp1*temp0 - DO n=1,numax - clff_u(n) = clff_u(n) + 2.0*gams_u(jc, n)*dyt/sref - cyff_u(n) = cyff_u(n) - 2.0*gams_u(jc, n)*dzt/sref - cdff_u(n) = cdff_u(n) + (gams_u(jc, n)*(dzt*vy-dyt*vz)+gams(jc - + )*(dzt*vy_u(n)-dyt*vz_u(n)))/sref - ENDDO - DO n=1,ncontrol - clff_d(n) = clff_d(n) + 2.0*gams_d(jc, n)*dyt/sref - cyff_d(n) = cyff_d(n) - 2.0*gams_d(jc, n)*dzt/sref - cdff_d(n) = cdff_d(n) + (gams_d(jc, n)*(dzt*vy-dyt*vz)+gams(jc - + )*(dzt*vy_d(n)-dyt*vz_d(n)))/sref - ENDDO - DO n=1,ndesign - clff_g(n) = clff_g(n) + 2.0*gams_g(jc, n)*dyt/sref - cyff_g(n) = cyff_g(n) - 2.0*gams_g(jc, n)*dzt/sref - cdff_g(n) = cdff_g(n) + (gams_g(jc, n)*(dzt*vy-dyt*vz)+gams(jc - + )*(dzt*vy_g(n)-dyt*vz_g(n)))/sref - ENDDO +C...Trefftz-plane drag is kinetic energy in crossflow +C + isurf = nsurfs(jc) + IF (lfload(isurf)) THEN +C-------add load from this strip only if it contributes to total load + temp1 = dyt/sref + clff_diff = clff_diff + 2.0*(temp1*gams_diff(jc)+gams(jc)*( + + dyt_diff-temp1*sref_diff)/sref) + clff = clff + 2.0*(gams(jc)*temp1) + temp1 = dzt/sref + cyff_diff = cyff_diff - 2.0*(temp1*gams_diff(jc)+gams(jc)*( + + dzt_diff-temp1*sref_diff)/sref) + cyff = cyff - 2.0*(gams(jc)*temp1) + temp0 = dzt*vy - dyt*vz + temp1 = gams(jc)/sref + cdff_diff = cdff_diff + temp0*(gams_diff(jc)-temp1*sref_diff)/ + + sref + temp1*(vy*dzt_diff+dzt*vy_diff-vz*dyt_diff-dyt* + + vz_diff) + cdff = cdff + temp1*temp0 + DO n=1,numax + clff_u(n) = clff_u(n) + 2.0*gams_u(jc, n)*dyt/sref + cyff_u(n) = cyff_u(n) - 2.0*gams_u(jc, n)*dzt/sref + cdff_u(n) = cdff_u(n) + (gams_u(jc, n)*(dzt*vy-dyt*vz)+gams( + + jc)*(dzt*vy_u(n)-dyt*vz_u(n)))/sref + ENDDO + DO n=1,ncontrol + clff_d(n) = clff_d(n) + 2.0*gams_d(jc, n)*dyt/sref + cyff_d(n) = cyff_d(n) - 2.0*gams_d(jc, n)*dzt/sref + cdff_d(n) = cdff_d(n) + (gams_d(jc, n)*(dzt*vy-dyt*vz)+gams( + + jc)*(dzt*vy_d(n)-dyt*vz_d(n)))/sref + ENDDO + DO n=1,ndesign + clff_g(n) = clff_g(n) + 2.0*gams_g(jc, n)*dyt/sref + cyff_g(n) = cyff_g(n) - 2.0*gams_g(jc, n)*dzt/sref + cdff_g(n) = cdff_g(n) + (gams_g(jc, n)*(dzt*vy-dyt*vz)+gams( + + jc)*(dzt*vy_g(n)-dyt*vz_g(n)))/sref + ENDDO + END IF ENDDO C C---- Double the X,Z forces, zero Y force for a Y symmetric case diff --git a/src/ad_src/reverse_ad_src/aero_b.f b/src/ad_src/reverse_ad_src/aero_b.f index a655209..28b94f8 100644 --- a/src/ad_src/reverse_ad_src/aero_b.f +++ b/src/ad_src/reverse_ad_src/aero_b.f @@ -1049,9 +1049,15 @@ SUBROUTINE SFFORC_B() C C******************************************************************* C -C...Store strip X,Y,Z body axes forces -C (these are normalized by strip area and moments are referred to -C c/4 point and are normalized by strip chord and area) +C--- At this point the forces are accumulated for the strip in body axes, +C referenced to strip 1/4 chord and normalized by strip area and chord +C CFX, CFY, CFZ ! body axes forces +C CMX, CMY, CMZ ! body axes moments about c/4 point +C CNC ! strip spanloading CN*chord +C CDV_LSTRP ! strip viscous drag in stability axes +C +C...Store strip X,Y,Z body axes forces and moments +C referred to c/4 and strip area and chord cf_strp(1, j) = cfx cf_strp(2, j) = cfy cf_strp(3, j) = cfz @@ -1059,12 +1065,14 @@ SUBROUTINE SFFORC_B() cm_strp(2, j) = cmy cm_strp(3, j) = cmz C -C...Transform strip body axes forces into stability axes - cdstrp(j) = cfx*cosa + cfz*sina - clstrp(j) = -(cfx*sina) + cfz*cosa +C...Strip body axes forces, referred to strip area and chord cxstrp(j) = cfx cystrp(j) = cfy czstrp(j) = cfz +C...Transform strip body axes forces into stability axes, +C referred to strip area and chord + cdstrp(j) = cfx*cosa + cfz*sina + clstrp(j) = -(cfx*sina) + cfz*cosa C cdst_a(j) = -(cfx*sina) + cfz*cosa clst_a(j) = -(cfx*cosa) - cfz*sina @@ -1093,11 +1101,12 @@ SUBROUTINE SFFORC_B() czst_g(j, n) = cfz_g(n) ENDDO C -C... Set strip moments about the overall moment reference point XYZREF -C (still normalized by strip area and chord) +C------ vector from chord c/4 reference point to case reference point XYZREF r(1) = xr - xyzref(1) r(2) = yr - xyzref(2) r(3) = zr - xyzref(3) +C... Strip moments in body axes about the case moment reference point XYZREF +C normalized by strip area and chord crstrp(j) = cmx + (cfz*r(2)-cfy*r(3))/cr cmstrp(j) = cmy + (cfx*r(3)-cfz*r(1))/cr cnstrp(j) = cmz + (cfy*r(1)-cfx*r(2))/cr @@ -1120,9 +1129,8 @@ SUBROUTINE SFFORC_B() cnst_g(j, n) = cmz_g(n) + (cfy_g(n)*r(1)-cfx_g(n)*r(2))/cr ENDDO C -C...Take components of X,Y,Z forces in local strip axes -C (axial/normal and lift/drag) -C in plane normal to (possibly dihedralled) strip +C...Components of X,Y,Z forces in local strip axes +C axial/normal forces and lift/drag in plane normal to dihedral of strip cl_lstrp(j) = ulift(1)*cfx + ulift(2)*cfy + ulift(3)*cfz cd_lstrp(j) = udrag(1)*cfx + udrag(2)*cfy + udrag(3)*cfz caxlstrp(j) = cfx @@ -1133,6 +1141,7 @@ SUBROUTINE SFFORC_B() rrot(1) = xsref(j) - xyzref(1) rrot(2) = ysref(j) - xyzref(2) rrot(3) = zsref(j) - xyzref(3) +C print *,"WROT ",WROT C C------ set total effective velocity = freestream + rotation CALL CROSS(rrot, wrot, vrot) @@ -1202,14 +1211,17 @@ SUBROUTINE SFFORC_B() cr = chord(j) C C -C -C +C--- Surface lift and drag referenced to case SREF, CREF, BREF +C--- Surface body axes forces referenced to case SREF, CREF, BREF +C--- Surface body axes moments referenced to case SREF, CREF, BREF about XYZREF C C--- Bug fix, HHY/S.Allmaras +C--- Surface viscous drag referenced to case SREF, CREF, BREF C ENDDO CALL PUSHINTEGER4(jj - 1) -C--- Surface hinge moments defined by surface LE moment about hinge vector +C +C--- Surface hinge moments defined by surface LE moment about hinge vector Ccc CMLE_SRF(IS) = DOT(CM_SRF(1,IS),VHINGE(1,IS)) C C @@ -1399,15 +1411,15 @@ SUBROUTINE SFFORC_B() cls_d_diff(is, n) = cls_d_diff(is, n) + cltot_d_diff(n) cds_d_diff(is, n) = cds_d_diff(is, n) + cdtot_d_diff(n) ENDDO - cdvsurf_diff(is) = cdvsurf_diff(is) + cdvtot_diff cnsurf_diff(is) = cnsurf_diff(is) + cntot_diff cmsurf_diff(is) = cmsurf_diff(is) + cmtot_diff crsurf_diff(is) = crsurf_diff(is) + crtot_diff + cdvsurf_diff(is) = cdvsurf_diff(is) + cdvtot_diff + clsurf_diff(is) = clsurf_diff(is) + cltot_diff + cdsurf_diff(is) = cdsurf_diff(is) + cdtot_diff czsurf_diff(is) = czsurf_diff(is) + cztot_diff cysurf_diff(is) = cysurf_diff(is) + cytot_diff cxsurf_diff(is) = cxsurf_diff(is) + cxtot_diff - clsurf_diff(is) = clsurf_diff(is) + cltot_diff - cdsurf_diff(is) = cdsurf_diff(is) + cdtot_diff END IF CALL POPINTEGER4(ad_to) DO jj=ad_to,1,-1 @@ -2180,18 +2192,18 @@ SUBROUTINE SFFORC_B() CALL PUSHCONTROL1B(1) END IF C -C... Set strip moments about the overall moment reference point XYZREF -C (still normalized by strip area and chord) +C------ vector from chord c/4 reference point to case reference point XYZREF CALL PUSHREAL8(r(1)) r(1) = xr - xyzref(1) CALL PUSHREAL8(r(2)) r(2) = yr - xyzref(2) CALL PUSHREAL8(r(3)) r(3) = zr - xyzref(3) +C... Strip moments in body axes about the case moment reference point XYZREF +C normalized by strip area and chord C -C...Take components of X,Y,Z forces in local strip axes -C (axial/normal and lift/drag) -C in plane normal to (possibly dihedralled) strip +C...Components of X,Y,Z forces in local strip axes +C axial/normal forces and lift/drag in plane normal to dihedral of strip C C------ vector at chord reference point from rotation axes CALL PUSHREAL8(rrot(1)) @@ -2200,6 +2212,7 @@ SUBROUTINE SFFORC_B() rrot(2) = ysref(j) - xyzref(2) CALL PUSHREAL8(rrot(3)) rrot(3) = zsref(j) - xyzref(3) +C print *,"WROT ",WROT C C------ set total effective velocity = freestream + rotation C @@ -2319,18 +2332,18 @@ SUBROUTINE SFFORC_B() clst_d_diff(j, n) = 0.D0 cdst_d_diff(j, n) = 0.D0 ENDDO - cfz_diff = cfz_diff + czstrp_diff(j) + cosa*clstrp_diff(j) + - + sina*cdstrp_diff(j) - czstrp_diff(j) = 0.D0 - cfy_diff = cfy_diff + cystrp_diff(j) - cystrp_diff(j) = 0.D0 - cfx_diff = cfx_diff + cxstrp_diff(j) + cosa*cdstrp_diff(j) - - + sina*clstrp_diff(j) - cxstrp_diff(j) = 0.D0 + cfz_diff = cfz_diff + cosa*clstrp_diff(j) + sina*cdstrp_diff(j) + + + czstrp_diff(j) cosa_diff = cosa_diff + cfz*clstrp_diff(j) + cfx*cdstrp_diff(j) + cfx_diff = cfx_diff + cosa*cdstrp_diff(j) - sina*clstrp_diff(j) + + + cxstrp_diff(j) sina_diff = sina_diff + cfz*cdstrp_diff(j) - cfx*clstrp_diff(j) clstrp_diff(j) = 0.D0 cdstrp_diff(j) = 0.D0 + czstrp_diff(j) = 0.D0 + cfy_diff = cfy_diff + cystrp_diff(j) + cystrp_diff(j) = 0.D0 + cxstrp_diff(j) = 0.D0 CALL POPCONTROL1B(branch) IF (branch .EQ. 0) THEN cdv_clv_diff = 0.D0 diff --git a/src/ad_src/reverse_ad_src/aic_b.f b/src/ad_src/reverse_ad_src/aic_b.f index c83eb8a..8745bf5 100644 --- a/src/ad_src/reverse_ad_src/aic_b.f +++ b/src/ad_src/reverse_ad_src/aic_b.f @@ -409,6 +409,8 @@ SUBROUTINE VORVELC_B(x, x_diff, y, y_diff, z, z_diff, lbound, x1, + rcore, rcore_diff) C---------------------------------------------------------- C Same as VORVEL, with finite core radius +C Uses Scully (also Burnham-Hallock) core model +C Vtan = Gam/2*pi . r/(r^2 +rcore^2) C---------------------------------------------------------- LOGICAL lbound C diff --git a/src/ad_src/reverse_ad_src/amake_b.f b/src/ad_src/reverse_ad_src/amake_b.f index 4364962..fe362e3 100644 --- a/src/ad_src/reverse_ad_src/amake_b.f +++ b/src/ad_src/reverse_ad_src/amake_b.f @@ -253,6 +253,7 @@ SUBROUTINE MAKESURF_B(isurf) REAL ypt1_diff REAL yscale REAL yscale_diff + INTEGER ii REAL width REAL width_diff REAL chordl @@ -332,6 +333,8 @@ SUBROUTINE MAKESURF_B(isurf) REAL fracte_diff INTRINSIC MAX INTRINSIC MIN + REAL zl + REAL zu REAL sum REAL wtot INTEGER jj @@ -401,7 +404,7 @@ SUBROUTINE MAKESURF_B(isurf) idx_strip = jfrst(isurf) C C----------------------------------------------------------------- -C---- section arc lengths of wing trace in y-z plane +C---- Arc length positions of sections in wing trace in y-z plane yzlen(1) = 0. DO isec=2,nsec(isurf) dy = xyzles(2, isec, isurf) - xyzles(2, isec-1, isurf) @@ -466,7 +469,13 @@ SUBROUTINE MAKESURF_B(isurf) END IF ELSE C -C----- set spanwise spacing using overall parameters NVS(ISURF), SSPACE +C +C----- Otherwise, set spanwise spacing using the SURFACE spanwise +C parameters NVS, SSPACE +C +C This spanwise spacing is modified (fudged) to align vortex edges +C with SECTIONs as defined. This allows CONTROLs to be defined +C without bridging vortex strips C nspace = 2*nvs(isurf) + 1 IF (nspace .GT. kpmax) THEN @@ -511,7 +520,7 @@ SUBROUTINE MAKESURF_B(isurf) CALL PUSHINTEGER4(isec) ad_count0 = 1 C -C----- fudge Glauert angles to make nodes match up exactly with interior sections +C----- fudge spacing array to make nodes match up exactly with interior sections DO isec=2,nsec(isurf)-1 CALL PUSHINTEGER4(ipt1) ipt1 = iptloc(isec-1) @@ -521,6 +530,7 @@ SUBROUTINE MAKESURF_B(isurf) GOTO 110 ELSE C +C----- fudge spacing to this section so that nodes match up exactly with section CALL PUSHREAL8(ypt1) ypt1 = ypt(ipt1) CALL PUSHREAL8(yscale) @@ -541,6 +551,7 @@ SUBROUTINE MAKESURF_B(isurf) CALL PUSHINTEGER4(ivs - 1) CALL PUSHINTEGER4(ad_from0) C +C----- check for unique spacing node for next section, if not we need more nodes CALL PUSHINTEGER4(ipt1) ipt1 = iptloc(isec) CALL PUSHINTEGER4(ipt2) @@ -549,6 +560,7 @@ SUBROUTINE MAKESURF_B(isurf) GOTO 120 ELSE C +C----- fudge spacing to this section so that nodes match up exactly with section CALL PUSHREAL8(ypt1) ypt1 = ypt(ipt1) CALL PUSHREAL8(yscale) @@ -584,7 +596,7 @@ SUBROUTINE MAKESURF_B(isurf) STOP END IF END IF -C +Cc#endif C C C==================================================== @@ -742,6 +754,9 @@ SUBROUTINE MAKESURF_B(isurf) END IF C C +Cc#ifdef USE_CPOML +C +Cc#endif CALL PUSHREAL8(chsin) chsin = chsinl + fc*(chsinr-chsinl) CALL PUSHREAL8(chcos) @@ -873,6 +888,8 @@ SUBROUTINE MAKESURF_B(isurf) C-------- go over vortices in this strip idx_vor = ijfrst(idx_strip) C NVOR = NVOR + 1 +C change all NVOR indices into idx_vor +C change all NSTRIP indices into idx_strip DO ivc=1,nvc(isurf) C C @@ -937,10 +954,19 @@ SUBROUTINE MAKESURF_B(isurf) ENDDO C C---------- TE control point used only if surface sheds a wake +C +Cc#ifdef USE_CPOML +C... nodal grid associated with vortex strip (aft-panel nodes) +C... NOTE: airfoil in plane of wing, but not rotated perpendicular to dihedral; +C... retained in (x,z) plane at this point +C +C +C +C +Cc#endif CALL PUSHINTEGER4(idx_vor) idx_vor = idx_vor + 1 ENDDO -C C CALL PUSHINTEGER4(idx_strip) idx_strip = idx_strip + 1 @@ -2167,6 +2193,9 @@ SUBROUTINE SDUPL_B(nn, ypt, msg) INTEGER klen INTRINSIC LEN INTEGER k + INTEGER isec + INTEGER idup + INTEGER iorg REAL yoff INTEGER idx_strip INTEGER ivs @@ -2287,6 +2316,9 @@ SUBROUTINE SDUPL_B(nn, ypt, msg) C C--- Image flag reversed (set to -IMAGS) for imaged surfaces C +Cc#ifdef USE_CPOML +Cc#endif +C C C--- Create image strips, to maintain the same sense of positive GAMMA C these have the 1 and 2 strip edges reversed (i.e. root is edge 2, @@ -2301,6 +2333,9 @@ SUBROUTINE SDUPL_B(nn, ypt, msg) jji = jfrst(nni) + ivs - 1 jj = jfrst(nn) + ivs - 1 C +Cc#ifdef USE_CPOML +C +Cc#endif C n = ndesign + 1 CALL PUSHINTEGER4(n - 1) @@ -2340,12 +2375,19 @@ SUBROUTINE SDUPL_B(nn, ypt, msg) rsgn = vrefl(jj, n) ENDDO CALL PUSHINTEGER4(n - 1) +C +Cc#ifdef USE_CPOML +C... nodal grid associated with vortex strip +C +Cc#endif idx_vor = idx_vor + 1 ad_count0 = ad_count0 + 1 END IF ENDDO CALL PUSHCONTROL1B(0) CALL PUSHINTEGER4(ad_count0) +C + idx_strip = idx_strip + 1 CALL PUSHINTEGER4(ivs) CALL PUSHCONTROL1B(1) END IF diff --git a/src/ad_src/reverse_ad_src/atpforc_b.f b/src/ad_src/reverse_ad_src/atpforc_b.f index c453fbe..2842965 100644 --- a/src/ad_src/reverse_ad_src/atpforc_b.f +++ b/src/ad_src/reverse_ad_src/atpforc_b.f @@ -79,6 +79,7 @@ SUBROUTINE TPFORC_B() REAL rsq1_diff REAL rsq2 REAL rsq2_diff + INTEGER isurf REAL ar REAL ar_diff REAL spanef_cl @@ -214,11 +215,19 @@ SUBROUTINE TPFORC_B() ENDDO C C +C C...Trefftz-plane drag is kinetic energy in crossflow C - clff = clff + 2.0*gams(jc)*dyt/sref - cyff = cyff - 2.0*gams(jc)*dzt/sref - cdff = cdff + gams(jc)*(dzt*vy-dyt*vz)/sref + isurf = nsurfs(jc) + IF (lfload(isurf)) THEN +C-------add load from this strip only if it contributes to total load + clff = clff + 2.0*gams(jc)*dyt/sref + cyff = cyff - 2.0*gams(jc)*dzt/sref + cdff = cdff + gams(jc)*(dzt*vy-dyt*vz)/sref + CALL PUSHCONTROL1B(0) + ELSE + CALL PUSHCONTROL1B(1) + END IF ENDDO C C---- Double the X,Z forces, zero Y force for a Y symmetric case @@ -273,22 +282,30 @@ SUBROUTINE TPFORC_B() ENDDO ENDDO DO jc=nstrip,1,-1 - dyt = rt2(2, jc) - rt1(2, jc) - dzt = rt2(3, jc) - rt1(3, jc) - temp = gams(jc)/sref - temp_diff1 = (dzt*vy-dyt*vz)*cdff_diff/sref - temp_diff0 = temp*cdff_diff - vy_diff = dzt*temp_diff0 - vz_diff = -(dyt*temp_diff0) - gams_diff(jc) = gams_diff(jc) + temp_diff1 + dyt*2.0*clff_diff/ - + sref - dzt*2.0*cyff_diff/sref - sref_diff = sref_diff - temp*temp_diff1 - temp_diff1 = -(gams(jc)*2.0*cyff_diff/sref) - dzt_diff = vy*temp_diff0 + temp_diff1 - sref_diff = sref_diff - dzt*temp_diff1/sref - temp_diff1 = gams(jc)*2.0*clff_diff/sref - dyt_diff = temp_diff1 - vz*temp_diff0 - sref_diff = sref_diff - dyt*temp_diff1/sref + CALL POPCONTROL1B(branch) + IF (branch .EQ. 0) THEN + dyt = rt2(2, jc) - rt1(2, jc) + dzt = rt2(3, jc) - rt1(3, jc) + temp = gams(jc)/sref + temp_diff1 = (dzt*vy-dyt*vz)*cdff_diff/sref + temp_diff0 = temp*cdff_diff + vy_diff = dzt*temp_diff0 + vz_diff = -(dyt*temp_diff0) + gams_diff(jc) = gams_diff(jc) + temp_diff1 + dyt*2.0*clff_diff + + /sref - dzt*2.0*cyff_diff/sref + sref_diff = sref_diff - temp*temp_diff1 + temp_diff1 = -(gams(jc)*2.0*cyff_diff/sref) + dzt_diff = vy*temp_diff0 + temp_diff1 + sref_diff = sref_diff - dzt*temp_diff1/sref + temp_diff1 = gams(jc)*2.0*clff_diff/sref + dyt_diff = temp_diff1 - vz*temp_diff0 + sref_diff = sref_diff - dyt*temp_diff1/sref + ELSE + dyt_diff = 0.D0 + vy_diff = 0.D0 + vz_diff = 0.D0 + dzt_diff = 0.D0 + END IF ycntr = rtc(2, jc) zcntr = rtc(3, jc) ycntr_diff = 0.D0 diff --git a/src/aero.f b/src/aero.f index d437ba7..66a15bc 100644 --- a/src/aero.f +++ b/src/aero.f @@ -822,9 +822,15 @@ SUBROUTINE SFFORC C C******************************************************************* C -C...Store strip X,Y,Z body axes forces -C (these are normalized by strip area and moments are referred to -C c/4 point and are normalized by strip chord and area) +C--- At this point the forces are accumulated for the strip in body axes, +C referenced to strip 1/4 chord and normalized by strip area and chord +C CFX, CFY, CFZ ! body axes forces +C CMX, CMY, CMZ ! body axes moments about c/4 point +C CNC ! strip spanloading CN*chord +C CDV_LSTRP ! strip viscous drag in stability axes +C +C...Store strip X,Y,Z body axes forces and moments +C referred to c/4 and strip area and chord CF_STRP(1,J) = CFX CF_STRP(2,J) = CFY CF_STRP(3,J) = CFZ @@ -832,12 +838,14 @@ SUBROUTINE SFFORC CM_STRP(2,J) = CMY CM_STRP(3,J) = CMZ C -C...Transform strip body axes forces into stability axes - CDSTRP(J) = CFX*COSA + CFZ*SINA - CLSTRP(J) = -CFX*SINA + CFZ*COSA +C...Strip body axes forces, referred to strip area and chord CXSTRP(J) = CFX CYSTRP(J) = CFY CZSTRP(J) = CFZ +C...Transform strip body axes forces into stability axes, +C referred to strip area and chord + CDSTRP(J) = CFX*COSA + CFZ*SINA + CLSTRP(J) = -CFX*SINA + CFZ*COSA C CDST_A(J) = -CFX*SINA + CFZ*COSA CLST_A(J) = -CFX*COSA - CFZ*SINA @@ -866,11 +874,12 @@ SUBROUTINE SFFORC CZST_G(J,N) = CFZ_G(N) END DO C -C... Set strip moments about the overall moment reference point XYZREF -C (still normalized by strip area and chord) +C------ vector from chord c/4 reference point to case reference point XYZREF R(1) = XR - XYZREF(1) R(2) = YR - XYZREF(2) R(3) = ZR - XYZREF(3) +C... Strip moments in body axes about the case moment reference point XYZREF +C normalized by strip area and chord CRSTRP(J) = CMX + (CFZ*R(2) - CFY*R(3))/CR CMSTRP(J) = CMY + (CFX*R(3) - CFZ*R(1))/CR CNSTRP(J) = CMZ + (CFY*R(1) - CFX*R(2))/CR @@ -893,9 +902,8 @@ SUBROUTINE SFFORC CNST_G(J,N) = CMZ_G(N) + (CFY_G(N)*R(1) - CFX_G(N)*R(2))/CR ENDDO C -C...Take components of X,Y,Z forces in local strip axes -C (axial/normal and lift/drag) -C in plane normal to (possibly dihedralled) strip +C...Components of X,Y,Z forces in local strip axes +C axial/normal forces and lift/drag in plane normal to dihedral of strip CL_LSTRP(J) = ULIFT(1)*CFX + ULIFT(2)*CFY + ULIFT(3)*CFZ CD_LSTRP(J) = UDRAG(1)*CFX + UDRAG(2)*CFY + UDRAG(3)*CFZ CAXLSTRP(J) = CFX @@ -906,6 +914,7 @@ SUBROUTINE SFFORC RROT(1) = XSREF(J) - XYZREF(1) RROT(2) = YSREF(J) - XYZREF(2) RROT(3) = ZSREF(J) - XYZREF(3) +c print *,"WROT ",WROT C C------ set total effective velocity = freestream + rotation CALL CROSS(RROT,WROT,VROT) @@ -1028,18 +1037,20 @@ SUBROUTINE SFFORC ENAVE(2) = ENAVE(2) + SR*ENSY(J) ENAVE(3) = ENAVE(3) + SR*ENSZ(J) C +C--- Surface lift and drag referenced to case SREF, CREF, BREF CDSURF(IS) = CDSURF(IS) + CDSTRP(J)*SR/SREF CLSURF(IS) = CLSURF(IS) + CLSTRP(J)*SR/SREF -C +C--- Surface body axes forces referenced to case SREF, CREF, BREF CXSURF(IS) = CXSURF(IS) + CXSTRP(J)*SR/SREF CYSURF(IS) = CYSURF(IS) + CYSTRP(J)*SR/SREF CZSURF(IS) = CZSURF(IS) + CZSTRP(J)*SR/SREF -C +C--- Surface body axes moments referenced to case SREF, CREF, BREF about XYZREF CRSURF(IS) = CRSURF(IS) + CRSTRP(J)*(SR/SREF)*(CR/BREF) CMSURF(IS) = CMSURF(IS) + CMSTRP(J)*(SR/SREF)*(CR/CREF) CNSURF(IS) = CNSURF(IS) + CNSTRP(J)*(SR/SREF)*(CR/BREF) C C--- Bug fix, HHY/S.Allmaras +C--- Surface viscous drag referenced to case SREF, CREF, BREF CDVSURF(IS) = CDVSURF(IS) + CDV_LSTRP(J)*(SR/SREF) C CDS_A(IS) = CDS_A(IS) + CDST_A(J)*SR/SREF @@ -1146,23 +1157,26 @@ SUBROUTINE SFFORC ENDIF CL_SRF(IS) = DOT(ULIFT,CF_SRF(1,IS)) CD_SRF(IS) = DOT(UDRAG,CF_SRF(1,IS)) -C--- Surface hinge moments defined by surface LE moment about hinge vector +C +C--- Surface hinge moments defined by surface LE moment about hinge vector ccc CMLE_SRF(IS) = DOT(CM_SRF(1,IS),VHINGE(1,IS)) C C C------------------------------------------------- IF(LFLOAD(IS)) THEN -C------- Total forces summed from surface forces... -C- normalized to configuration reference quantities - CDTOT = CDTOT + CDSURF(IS) - CLTOT = CLTOT + CLSURF(IS) +C--- Total forces summed from surface forces +C normalized to case reference quantities SREF, CREF, BREF CXTOT = CXTOT + CXSURF(IS) CYTOT = CYTOT + CYSURF(IS) CZTOT = CZTOT + CZSURF(IS) + CDTOT = CDTOT + CDSURF(IS) + CLTOT = CLTOT + CLSURF(IS) + CDVTOT = CDVTOT + CDVSURF(IS) +C--- Total body axes moments about XYZREF summed from surface moments +C normalized to case reference quantities SREF, CREF, BREF CRTOT = CRTOT + CRSURF(IS) CMTOT = CMTOT + CMSURF(IS) CNTOT = CNTOT + CNSURF(IS) - CDVTOT = CDVTOT + CDVSURF(IS) C CDTOT_A = CDTOT_A + CDS_A(IS) CLTOT_A = CLTOT_A + CLS_A(IS) @@ -1456,7 +1470,6 @@ SUBROUTINE VINFAB C INCLUDE 'AVL.INC' C - SINA = SIN(ALFA) COSA = COS(ALFA) SINB = SIN(BETA) diff --git a/src/aic.f b/src/aic.f index 214e1ab..ee4656a 100644 --- a/src/aic.f +++ b/src/aic.f @@ -362,7 +362,8 @@ SUBROUTINE SRDSET(BETM,XYZREF, & SRC_U,DBL_U ) C---------------------------------------------------------- C Sets strengths of source+doublet line segments -C for each unit freestream and rotation component +C for 6 "unit" flow components consisting of: +C 3 unit (X,Y,Z) freestream and 3 unit (X,Y,Z) rotations C C Input C ----- @@ -567,6 +568,8 @@ SUBROUTINE VORVELC(X,Y,Z,LBOUND,X1,Y1,Z1,X2,Y2,Z2,BETA, & U,V,W, RCORE) C---------------------------------------------------------- C Same as VORVEL, with finite core radius +C Uses Scully (also Burnham-Hallock) core model +C Vtan = Gam/2*pi . r/(r^2 +rcore^2) C---------------------------------------------------------- LOGICAL LBOUND C diff --git a/src/ainput.f b/src/ainput.f index 5e9792f..8911cc2 100644 --- a/src/ainput.f +++ b/src/ainput.f @@ -1,8 +1,8 @@ C*********************************************************************** C Module: ainput.f -C +C C Copyright (C) 2002 Mark Drela, Harold Youngren -C +C C This program is free software; you can redistribute it and/or modify C it under the terms of the GNU General Public License as published by C the Free Software Foundation; either version 2 of the License, or @@ -28,8 +28,8 @@ SUBROUTINE INPUT(LUN,FNAME,FERR) LOGICAL FERR C CHARACTER*4 KEYWD - CHARACTER*120 CNAME, ANAME - CHARACTER*128 LINE + CHARACTER*256 CNAME, ANAME + CHARACTER*256 LINE LOGICAL LHINGE C REAL CLX(3), CDX(3) @@ -40,7 +40,6 @@ SUBROUTINE INPUT(LUN,FNAME,FERR) REAL XB(IBX), YB(IBX) REAL XIN(IBX), YIN(IBX), TIN(IBX) REAL XBOD(IBX), YBOD(IBX), TBOD(IBX) - C C---- max number of control or design variable declaration lines per section C @@ -53,13 +52,13 @@ SUBROUTINE INPUT(LUN,FNAME,FERR) C C---------------------------------------------------- FERR = .FALSE. -C +C OPEN(UNIT=LUN,FILE=FNAME,STATUS='OLD',ERR=3) GO TO 6 C 3 CONTINUE CALL STRIP(FNAME,NFN) - WRITE(*,*) + WRITE(*,*) WRITE(*,*) '** Open error on file: ', FNAME(1:NFN) FNAME = FNAME(1:NFN) // '.avl' WRITE(*,*) ' Trying alternative: ', FNAME(1:NFN+4) @@ -68,7 +67,7 @@ SUBROUTINE INPUT(LUN,FNAME,FERR) C 4 CONTINUE CALL STRIP(FNAME,NFN) - WRITE(*,*) + WRITE(*,*) WRITE(*,*) '** Open error on file: ', FNAME(1:NFN) FERR = .TRUE. RETURN @@ -99,6 +98,12 @@ SUBROUTINE INPUT(LUN,FNAME,FERR) C NCONTROL = 0 NDESIGN = 0 +C---- clear surface polar data + DO K=1,NFMAX + DO L=1, 6 + CLCDSRF(L,k) = 0.0 + END DO + END DO C LVISC = .FALSE. C @@ -361,20 +366,29 @@ SUBROUTINE INPUT(LUN,FNAME,FERR) ENDIF C CALL RDLINE(LUN,LINE,NLINE,ILINE) - C IF (ISURF.NE.0) THEN CC WRITE(*,*) ' + duplicate surface ',STITLE(ISURF) LDUPL(ISURF) = .TRUE. READ (LINE,*,ERR=990) YDUPL(ISURF) + IF(IYSYM.NE.0 .AND. YDUPL(ISURF).EQ.0.0) THEN + WRITE(*,*) 'ERROR: Redundant y-symmetry specifications...' + WRITE(*,*) ' IYSYM /= 0' + WRITE(*,*) ' YDUPLICATE 0.0' + WRITE(*,*) 'Can use one or the other, but not both.' + STOP + ENDIF ELSEIF(IBODY.NE.0) THEN CC WRITE(*,*) ' + duplicate body ',BTITLE(IBODY) LDUPL_B(IBODY) = .TRUE. READ (LINE,*,ERR=990) YDUPL_B(IBODY) - ENDIF -C - IF(IYSYM.NE.0) THEN - WRITE(*,*)'** Warning: Y-duplicate AND Y-sym specified' + IF(IYSYM.NE.0 .AND. YDUPL_B(IBODY).EQ.0.0) THEN + WRITE(*,*) 'ERROR: Redundant y-symmetry specifications...' + WRITE(*,*) ' IYSYM /= 0' + WRITE(*,*) ' YDUPLICATE 0.0' + WRITE(*,*) 'Can use one or the other, but not both.' + STOP + ENDIF ENDIF C C=========================================================================== @@ -445,6 +459,8 @@ SUBROUTINE INPUT(LUN,FNAME,FERR) C CALL RDLINE(LUN,LINE,NLINE,ILINE) READ (LINE,*,ERR=990) ADDINC(ISURF) + ! Don't do the conversion inplace! modified the other routines to + ! include conversion ! ADDINC(ISURF) = ADDINC(ISURF)*DTR C @@ -555,9 +571,20 @@ SUBROUTINE INPUT(LUN,FNAME,FERR) SASEC(2,ISEC, ISURF) = 0.0 TASEC(1,ISEC, ISURF) = 0.0 TASEC(2,ISEC, ISURF) = 0.0 -C ...no polar data +cc#ifdef USE_CPOML +C + XLASEC(1,ISEC,ISURF) = 0 + XLASEC(2,ISEC,ISURF) = 0 + XUASEC(1,ISEC,ISURF) = 0 + XUASEC(2,ISEC,ISURF) = 0 + ZLASEC(1,ISEC,ISURF) = 0 + ZLASEC(2,ISEC,ISURF) = 0 + ZUASEC(1,ISEC,ISURF) = 0 + ZUASEC(2,ISEC,ISURF) = 0 +cc#endif +C ...copy polar data from surface polar (if specified) DO L=1, 6 - CLCDSEC(L,ISEC,ISURF) = 0.0 + CLCDSEC(L,ISEC,ISURF) = CLCDSRF(L,ISURF) END DO C ...no control NSCON(ISEC, ISURF) = 0 @@ -568,7 +595,7 @@ SUBROUTINE INPUT(LUN,FNAME,FERR) CLAF(ISEC,ISURF) = 1.0 C C=========================================================================== - ELSEIF (KEYWD.EQ.'NACA') THEN + ELSEIF (KEYWD.EQ.'NACA') THEN C------ input NACA camberline C IF(ISURF.EQ.0 .OR. ISEC.EQ.0) THEN @@ -577,18 +604,35 @@ SUBROUTINE INPUT(LUN,FNAME,FERR) GO TO 10 ENDIF C - CALL RDLINE(LUN,LINE,NLINE,ILINE) -C - IB = INDEX(LINE,' ') - READ(LINE(1:IB-1),*,ERR=990) IDES - IF(LINE(IB:NLINE).NE.' ') THEN - READ(LINE(IB:NLINE),*,ERR=990) XFMIN, XFMAX -ccc WRITE(*,*) ' Using data in normalized range ',XFMIN,XFMAX - ELSE +C------ find blank (if any) after keyword + IB = INDEX(LINE(1:NLINE),' ') + + IF(IB.EQ.0 .OR. IB.GE.NLINE-2) THEN +C------- default airfoil x/c limits are the whole airfoil XFMIN = 0. XFMAX = 1. + ELSE +C------- read optional airfoil x/c limits + NINPUT = 2 + CALL GETFLT(LINE(IB:NLINE),RINPUT,NINPUT,ERROR) + IF(ERROR .OR. NINPUT.LT.2) THEN + XFMIN = 0. + XFMAX = 1. + ELSE + XFMIN = MAX( 0.0 , RINPUT(1) ) + XFMAX = MIN( 1.0 , RINPUT(2) ) + ENDIF ENDIF C + IF((XFMIN .GT. 0.01) .OR. (XFMAX .LT. 0.99)) THEN + LRANGE(ISURF) = .FALSE. + ENDIF +C +C------ read NACA 4-digit designator + CALL RDLINE(LUN,LINE,NLINE,ILINE) + READ(LINE,*,ERR=990) IDES + +C------ parse designator and set camber coordinates ICAM = IDES/1000 IPOS = (IDES-1000*ICAM)/100 ITHK = IDES-1000*ICAM-100*IPOS @@ -601,26 +645,43 @@ SUBROUTINE INPUT(LUN,FNAME,FERR) XF = XFMIN + & (XFMAX-XFMIN)*FLOAT(I-1)/FLOAT(NASEC(ISEC,ISURF)-1) - XASEC(I,ISEC, ISURF) = XF - TASEC(I,ISEC, ISURF) = - & ( 0.29690*SQRT(XF) + IF(XF .LT. P) THEN + SLP = C * 2.0*(P - XF) / P**2 + ELSEIF(XF .GT. P) THEN + SLP = C * 2.0*(P - XF) / (1.0-P)**2 + ELSE + SLP = 0 + ENDIF + THK = (0.29690*SQRT(XF) & - 0.12600*XF & - 0.35160*XF**2 & + 0.28430*XF**3 - & - 0.10150*XF**4 ) * T * 10.0 - IF (XF .LT. P) THEN - SASEC(I,ISEC,ISURF) = C * 2.0*(P - XF) / P**2 - ELSEIF(XF .GE. P) THEN - SASEC(I,ISEC,ISURF) = C * 2.0*(P - XF) / (1.0-P)**2 + & - 0.10150*XF**4) * T * 10.0 +C + XASEC(I,ISEC,ISURF) = XF + SASEC(I,ISEC,ISURF) = SLP + TASEC(I,ISEC,ISURF) = THK +cc#ifdef USE_CPOML +c... lower/upper coordinates (for oml) + IF(XF .LT. P) THEN + ZF = C * (2*P*XF - 1)*XF / P**2 + ELSEIF(XF .GT. P) THEN + ZF = C * ((1 - 2*P) + (2*P - XF)*XF) / (1-P)**2 ELSE - SASEC(I,ISEC,ISURF) = 0. + ZF = 0 ENDIF + TH = ATAN(SLP) + XLASEC(I,ISEC,ISURF) = XF + 0.5*THK*SIN(TH) + XUASEC(I,ISEC,ISURF) = XF - 0.5*THK*SIN(TH) + ZLASEC(I,ISEC,ISURF) = ZF - 0.5*THK*COS(TH) + ZUASEC(I,ISEC,ISURF) = ZF + 0.5*THK*COS(TH) +cc#endif ENDDO CALL NRMLIZ(NASEC(ISEC,ISURF),XASEC(1,ISEC,ISURF)) C C=========================================================================== - ELSE IF (KEYWD.EQ.'AIRF') THEN + ELSE IF (KEYWD.EQ.'AIRF') THEN C------ input y(x) for an airfoil, get camber then slopes via spline C IF(ISURF.EQ.0 .OR. ISEC.EQ.0) THEN @@ -628,12 +689,32 @@ SUBROUTINE INPUT(LUN,FNAME,FERR) WRITE(*,*) '** No active section for this airfoil' GO TO 10 ENDIF + +C------ find blank (if any) after keyword + IB = INDEX(LINE(1:NLINE),' ') + + IF(IB.EQ.0 .OR. IB.GE.NLINE-2) THEN +C------- default airfoil x/c limits are the whole airfoil + XFMIN = 0. + XFMAX = 1. + ELSE +C------- read optional airfoil x/c limits + NINPUT = 2 + CALL GETFLT(LINE(IB:NLINE),RINPUT,NINPUT,ERROR) + IF(ERROR .OR. NINPUT.LT.2) THEN + XFMIN = 0. + XFMAX = 1. + ELSE + XFMIN = MAX( 0.0 , RINPUT(1) ) + XFMAX = MIN( 1.0 , RINPUT(2) ) + ENDIF + ENDIF C - CALL RDLINE(LUN,LINE,NLINE,ILINE) -C - IB = INDEX(LINE,' ') - READ(LINE(IB:NLINE),*,ERR=990) XFMIN, XFMAX + IF((XFMIN .GT. 0.01) .OR. (XFMAX .LT. 0.99)) THEN + LRANGE(ISURF) = .FALSE. + ENDIF C +C------ read airfoil coordinates DO I = 1, 123456 IB = MIN(I,IBX) C @@ -651,7 +732,7 @@ SUBROUTINE INPUT(LUN,FNAME,FERR) C 40 CONTINUE IF(I.GT.IBX) THEN - WRITE(*,*) + WRITE(*,*) & '*** AINPUT: Airfoil array overflow. Increase IBX to', I STOP ENDIF @@ -670,6 +751,14 @@ SUBROUTINE INPUT(LUN,FNAME,FERR) & SASEC(I,ISEC,ISURF)) CALL AKIMA(XIN,TIN,NIN,XASEC(I,ISEC,ISURF), & TASEC(I,ISEC,ISURF),DUMMY) +cc#ifdef USE_CPOML +C... lower/upper coordinates (for oml) + CALL AKIMA(XIN,YIN,NIN,XASEC(I,ISEC,ISURF),ZC,DUMMY) + XLASEC(I,ISEC,ISURF) = XASEC(I,ISEC,ISURF) + XUASEC(I,ISEC,ISURF) = XASEC(I,ISEC,ISURF) + ZLASEC(I,ISEC,ISURF) = ZC - 0.5*TASEC(I,ISEC,ISURF) + ZUASEC(I,ISEC,ISURF) = ZC + 0.5*TASEC(I,ISEC,ISURF) +cc#endif END DO CALL NRMLIZ(NASEC(ISEC,ISURF),XASEC(1,ISEC,ISURF)) C @@ -677,7 +766,7 @@ SUBROUTINE INPUT(LUN,FNAME,FERR) GO TO 11 C C=========================================================================== - ELSEIF (KEYWD.EQ.'AFIL') THEN + ELSEIF (KEYWD.EQ.'AFIL') THEN C------ input y(x) from an airfoil coordinate file C IF(ISURF.EQ.0 .OR. ISEC.EQ.0) THEN @@ -685,33 +774,50 @@ SUBROUTINE INPUT(LUN,FNAME,FERR) WRITE(*,*) '** No active section for this airfoil' GO TO 10 ENDIF + +C------ find blank (if any) after keyword + IB = INDEX(LINE(1:NLINE),' ') + + IF(IB.EQ.0 .OR. IB.GE.NLINE-2) THEN +C------- default airfoil x/c limits are the whole airfoil + XFMIN = 0. + XFMAX = 1. + ELSE +C------- read optional airfoil x/c limits + NINPUT = 2 + CALL GETFLT(LINE(IB:NLINE),RINPUT,NINPUT,ERROR) + IF(ERROR .OR. NINPUT.LT.2) THEN + XFMIN = 0. + XFMAX = 1. + ELSE + XFMIN = MAX( 0.0 , RINPUT(1) ) + XFMAX = MIN( 1.0 , RINPUT(2) ) + ENDIF + ENDIF C + IF((XFMIN .GT. 0.01) .OR. (XFMAX .LT. 0.99)) THEN + LRANGE(ISURF) = .FALSE. + ENDIF +C +C------ read filename CALL RDLINE(LUN,LINE,NLINE,ILINE) -C---- parse file name and optional parameters -C double quotes checked to delimit file name to allow blanks in name - IDQ1 = INDEX(LINE,'"') + +C------ parse file name +C double quotes checked to delimit file name to allow blanks in name + IDQ1 = INDEX(LINE(1:NLINE),'"') IF(IDQ1.NE.0) THEN - IDQ2 = INDEX(LINE(IDQ1+1:),'"') - IF(IDQ2.GT.1) THEN - CNAME = LINE(IDQ1+1:IDQ2+IDQ1-1) - IB = IDQ2 + IDQ1 + 1 + IDQ2 = INDEX(LINE(IDQ1+1:NLINE),'"') + IF(IDQ2.GT.1) THEN + CNAME = LINE(IDQ1+1:IDQ2+IDQ1-1) ELSE - WRITE(*,9000) '** Bad quotes in file name ',ILINE, - & LINE(1:NLINE) - GO TO 10 + WRITE(*,9000) '** Bad quotes in file name ',ILINE, + & LINE(1:NLINE) + GO TO 10 ENDIF + ELSE -C---- Find blank after filename as delimiter for optional parameters - IB = INDEX(LINE,' ') - CNAME = LINE(1:IB) - ENDIF -C - IF(LINE(IB:NLINE).NE.' ') THEN - READ(LINE(IB:NLINE),*,ERR=990) XFMIN, XFMAX -CCC WRITE(*,*) ' Using data in normalized range ',XFMIN,XFMAX - ELSE - XFMIN = 0. - XFMAX = 1. + CNAME = LINE(1:NLINE) + ENDIF C CALL STRIP(CNAME,NCN) @@ -728,13 +834,18 @@ SUBROUTINE INPUT(LUN,FNAME,FERR) WRITE(*,*) '** Airfoil file not found : ',CNAME(1:NCN) WRITE(*,*) '** Using default zero-camber airfoil' end if - C C NASEC(ISEC,ISURF) = MIN( 50 , IBX ) DO I = 1, NASEC(ISEC,ISURF) XASEC(I,ISEC,ISURF) = FLOAT(I-1)/FLOAT(NASEC(ISEC,ISURF)-1) SASEC(I,ISEC,ISURF) = 0.0 TASEC(I,ISEC,ISURF) = 0.0 +cc#ifdef USE_CPOML + XLASEC(I,ISEC,ISURF) = 0 + XUASEC(I,ISEC,ISURF) = 0 + ZLASEC(I,ISEC,ISURF) = 0 + ZUASEC(I,ISEC,ISURF) = 0 +cc#endif ENDDO C ELSE @@ -752,12 +863,20 @@ SUBROUTINE INPUT(LUN,FNAME,FERR) & MMY,SASEC(I,ISEC,ISURF)) CALL AKIMA(XIN,TIN,NIN,XASEC(I,ISEC,ISURF), & TASEC(I,ISEC,ISURF),DUMMY) +cc#ifdef USE_CPOML +C... lower/upper coordinates (for oml) + CALL AKIMA(XIN,YIN,NIN,XASEC(I,ISEC,ISURF),ZC,DUMMY) + XLASEC(I,ISEC,ISURF) = XASEC(I,ISEC,ISURF) + XUASEC(I,ISEC,ISURF) = XASEC(I,ISEC,ISURF) + ZLASEC(I,ISEC,ISURF) = ZC - 0.5*TASEC(I,ISEC,ISURF) + ZUASEC(I,ISEC,ISURF) = ZC + 0.5*TASEC(I,ISEC,ISURF) +cc#endif END DO CALL NRMLIZ (NASEC(ISEC,ISURF),XASEC(1,ISEC,ISURF)) ENDIF C C=========================================================================== - ELSEIF (KEYWD.EQ.'BFIL') THEN + ELSEIF (KEYWD.EQ.'BFIL') THEN C------ input r(x) from an airfoil coordinate file C IF(IBODY.EQ.0) THEN @@ -765,14 +884,39 @@ SUBROUTINE INPUT(LUN,FNAME,FERR) WRITE(*,*) '** No active body for this shape' GO TO 10 ENDIF +C +C------ find blank (if any) after keyword + IB = INDEX(LINE(1:NLINE),' ') + + IF(IB.EQ.0 .OR. IB.GE.NLINE-2) THEN +C------- default x/c limits are the whole airfoil + XFMIN = 0. + XFMAX = 1. + ELSE +C------- read optional x/c limits + NINPUT = 2 + CALL GETFLT(LINE(IB:NLINE),RINPUT,NINPUT,ERROR) + IF(ERROR .OR. NINPUT.LT.2) THEN + XFMIN = 0. + XFMAX = 1. + ELSE + XFMIN = MAX( 0.0 , RINPUT(1) ) + XFMAX = MIN( 1.0 , RINPUT(2) ) + ENDIF + ENDIF +C + IF((XFMIN .GT. 0.01) .OR. (XFMAX .LT. 0.99)) THEN + LRANGE(ISURF) = .FALSE. + ENDIF C CALL RDLINE(LUN,LINE,NLINE,ILINE) + C---- parse file name and optional parameters -C double quotes checked to delimit file name to allow blanks in name - IDQ1 = INDEX(LINE,'"') +C double quotes checked to delimit file name to allow blanks in name + IDQ1 = INDEX(LINE(1:NLINE),'"') IF(IDQ1.NE.0) THEN - IDQ2 = INDEX(LINE(IDQ1+1:),'"') - IF(IDQ2.GT.1) THEN + IDQ2 = INDEX(LINE(IDQ1+1:NLINE),'"') + IF(IDQ2.GT.1) THEN CNAME = LINE(IDQ1+1:IDQ2+IDQ1-1) IB = IDQ2 + IDQ1 + 1 ELSE @@ -781,19 +925,11 @@ SUBROUTINE INPUT(LUN,FNAME,FERR) GO TO 10 ENDIF ELSE -C---- Find blank after filename as delimiter for optional parameters - IB = INDEX(LINE,' ') - CNAME = LINE(1:IB) - ENDIF -C - IF(LINE(IB:NLINE).NE.' ') THEN - READ(LINE(IB:NLINE),*,ERR=990) XFMIN, XFMAX -CCC WRITE(*,*) ' Using data in normalized range ',XFMIN,XFMAX - ELSE - XFMIN = 0. - XFMAX = 1. + CNAME = LINE(1:NLINE) + ENDIF -C + + CALL STRIP(CNAME,NCN) if (lverbose) then WRITE(*,*) ' Reading body shape from file: ',CNAME(1:NCN) @@ -803,58 +939,105 @@ SUBROUTINE INPUT(LUN,FNAME,FERR) CALL READBL(CNAME,IBX,NBLDS,XB,YB,NB,NBL, & ANAME,XINL,XOUT,YBOT,YTOP) C + IF(NB.LE.2) THEN + WRITE(*,*) '** Error reading body from ', CNAME(1:NCN) + WRITE(*,*) ' Body not defined' + IBODY = 0 + NBODY = NBODY - 1 + GO TO 10 + ELSE C------ set thread line y, and thickness t ( = 2r) - NBOD = MIN( 50 , IBX ) - CALL GETCAM(XB,YB,NB,XBOD,YBOD,TBOD,NBOD,.FALSE.) + NBOD = MIN( 50 , IBX ) + CALL GETCAM(XB,YB,NB,XBOD,YBOD,TBOD,NBOD,.FALSE.) + ENDIF C C=========================================================================== - ELSEIF (KEYWD.EQ.'CDCL') THEN + ELSEIF (KEYWD.EQ.'CDCL') THEN C------ input approximate CD(CL) polar defining data C - IF(ISURF.EQ.0 .OR. ISEC.EQ.0) THEN + IF(ISURF.EQ.0 .AND. ISEC.EQ.0) THEN WRITE(*,9000) '** Misplaced line', ILINE, LINE(1:NLINE) - WRITE(*,*) '** No active section for this polar' + WRITE(*,*) '** No active surface or section for this polar' GO TO 10 ENDIF - +C +C------ if defining surface before sections store polar for surface + IF(ISURF.NE.0 .AND. ISEC.EQ.0) THEN + if (lverbose) then + WRITE(*,*) '** Polar data for surface (all sections)' + ENDIF CALL RDLINE(LUN,LINE,NLINE,ILINE) READ(LINE,*,ERR=990) CLX(1),CDX(1),CLX(2),CDX(2),CLX(3),CDX(3) C LMAX = 1 LMIN = 1 - DO L = 2, 3 + DO L = 2, 3 IF(CLX(L).GT.CLX(LMAX)) LMAX = L IF(CLX(L).LT.CLX(LMIN)) LMIN = L END DO C - IF(ISEC.GT.1) THEN - IF(CLCDSEC(4,ISEC-1,ISURF).LE.0.0) THEN - WRITE(*,*) '* AINPUT: previous section defined with no polar' +C------ Trick: sum must be 6 so we can get the "other" index + LMID = 6 - (LMIN+LMAX) + CLCDSRF(1,ISURF) = CLX(LMIN) + CLCDSRF(2,ISURF) = CDX(LMIN) + CLCDSRF(3,ISURF) = CLX(LMID) + CLCDSRF(4,ISURF) = CDX(LMID) + CLCDSRF(5,ISURF) = CLX(LMAX) + CLCDSRF(6,ISURF) = CDX(LMAX) + if (lverbose) then + WRITE(*,1700) CLX(LMIN),CDX(LMIN), + & CLX(LMID),CDX(LMID), + & CLX(LMAX),CDX(LMAX) + endif + 1700 FORMAT(' Reading CD(CL) data for surface', + & /' CLneg = ',F8.3,' CD@CLneg = ',F10.5, + & /' CL@CDmin = ',F8.3,' CDmin = ',F10.5, + & /' CLpos = ',F8.3,' CD@CLpos = ',F10.5) + LVISC = .TRUE. +C +C------ define polar for the current active section + ELSEIF(ISURF.NE.0 .AND. ISEC.NE.0) THEN +C + CALL RDLINE(LUN,LINE,NLINE,ILINE) + READ(LINE,*,ERR=990) CLX(1),CDX(1),CLX(2),CDX(2),CLX(3),CDX(3) +C + LMAX = 1 + LMIN = 1 + DO L = 2, 3 + IF(CLX(L).GT.CLX(LMAX)) LMAX = L + IF(CLX(L).LT.CLX(LMIN)) LMIN = L + END DO +C + IF(ISEC.GT.1) THEN + IF(CLCDSEC(4,ISEC-1,ISURF).LE.0.0) THEN + WRITE(*,*) '* AINPUT: previous section defined with no polar' + ENDIF ENDIF - ENDIF C C------ Trick: sum must be 6 so we can get the "other" index LMID = 6 - (LMIN+LMAX) - CLCDSEC(1,ISEC,ISUF) = CLX(LMIN) - CLCDSEC(2,ISEC,ISUF) = CDX(LMIN) - CLCDSEC(3,ISEC,ISUF) = CLX(LMID) - CLCDSEC(4,ISEC,ISUF) = CDX(LMID) - CLCDSEC(5,ISEC,ISUF) = CLX(LMAX) - CLCDSEC(6,ISEC,ISUF) = CDX(LMAX) + CLCDSEC(1,ISEC,ISURF) = CLX(LMIN) + CLCDSEC(2,ISEC,ISURF) = CDX(LMIN) + CLCDSEC(3,ISEC,ISURF) = CLX(LMID) + CLCDSEC(4,ISEC,ISURF) = CDX(LMID) + CLCDSEC(5,ISEC,ISURF) = CLX(LMAX) + CLCDSEC(6,ISEC,ISURF) = CDX(LMAX) if (LVERBOSE) then - WRITE(*,1700) CLX(LMIN),CDX(LMIN), + WRITE(*,1701) CLX(LMIN),CDX(LMIN), & CLX(LMID),CDX(LMID), & CLX(LMAX),CDX(LMAX) -1700 FORMAT(' Reading CD(CL) data for section', + endif + 1701 FORMAT(' Reading CD(CL) data for section', & /' CLneg = ',F8.3,' CD@CLneg = ',F10.5, & /' CL@CDmin = ',F8.3,' CDmin = ',F10.5, & /' CLpos = ',F8.3,' CD@CLpos = ',F10.5) + LVISC = .TRUE. +C ENDIF - LVISC = .TRUE. C C=========================================================================== - ELSEIF (KEYWD.EQ.'CLAF') THEN + ELSEIF (KEYWD.EQ.'CLAF') THEN C------ input dCL/da scaling factor C IF(ISURF.EQ.0 .OR. ISEC.EQ.0) THEN @@ -899,7 +1082,10 @@ SUBROUTINE INPUT(LUN,FNAME,FERR) C C------ see if this control variable has already been declared DO N = 1, NCONTROL - IF(LINE(1:NNAME) .EQ. DNAME(N)(1:NNAME)) THEN + NDNAME = INDEX(DNAME(N),' ') - 1 ! added 18 Nov 2021 MD + IF(NNAME .EQ. NDNAME .AND. ! added 18 Nov 2021 MD + & LINE(1:NNAME) .EQ. DNAME(N)(1:NNAME)) THEN ! modified 18 Nov 2021 +ccc IF(LINE(1:NNAME) .EQ. DNAME(N)(1:NNAME)) THEN ICONTROL = N GO TO 62 ENDIF @@ -950,7 +1136,7 @@ SUBROUTINE INPUT(LUN,FNAME,FERR) ENDIF C C=========================================================================== - ELSEIF (KEYWD.EQ.'DESI') THEN + ELSEIF (KEYWD.EQ.'DESI') THEN C------ link section to design variable and weight C IF(ISURF.EQ.0 .OR. ISEC.EQ.0) THEN @@ -968,7 +1154,7 @@ SUBROUTINE INPUT(LUN,FNAME,FERR) C------ extract design name NNAME = INDEX(LINE,' ') - 1 IF(NNAME.LE.0) THEN - WRITE(*,9000) ' *** Bad design declaration line', + WRITE(*,9000) ' *** Bad design declaration line', & ILINE, LINE(1:NLINE) STOP ENDIF @@ -1015,7 +1201,10 @@ SUBROUTINE INPUT(LUN,FNAME,FERR) C C********************************************************************* C -C WRITE (*,2018) MACH0,NBODY,NSURF,NSTRIP,NVOR + if (lverbose) then + WRITE (*,2018) MACH0,NBODY,NSURF,NSTRIP,NVOR + WRITE (*,2019) NCONTROL,NDESIGN + end if C IF(IYSYM.GT.0) WRITE (*,2024) YSYM IF(IYSYM.LT.0) WRITE (*,2025) YSYM @@ -1023,10 +1212,13 @@ SUBROUTINE INPUT(LUN,FNAME,FERR) IF(IZSYM.LT.0) WRITE (*,2027) ZSYM C 2018 FORMAT (/' Mach =',F10.4,' (default)' - & /' Nbody =',I4,5X, - & ' Nsurf =',I4,5X, - & ' Nstrp =',I4,5X, - & ' Nvor =',I4) + & //1X,I4,' Bodies' + & /1X,I4,' Solid surfaces' + & /1X,I4,' Strips' + & /1X,I4,' Vortices') + 2019 FORMAT (/1X,I4,' Control variables' + & /1X,I4,' Design parameters') + 2024 FORMAT (/' Y Symmetry: Wall plane at Ysym =',F10.4) 2025 FORMAT (/' Y Symmetry: Free surface at Ysym =',F10.4) 2026 FORMAT (/' Z Symmetry: Ground plane at Zsym =',F10.4) @@ -1101,6 +1293,6 @@ SUBROUTINE TOUPER(INPUT) 10 CONTINUE C RETURN - END + END diff --git a/src/amake.f b/src/amake.f index 29ea45c..9d9ff48 100644 --- a/src/amake.f +++ b/src/amake.f @@ -92,7 +92,7 @@ SUBROUTINE MAKESURF(ISURF) idx_strip = JFRST(ISURF) C C----------------------------------------------------------------- -C---- section arc lengths of wing trace in y-z plane +C---- Arc length positions of sections in wing trace in y-z plane YZLEN(1) = 0. DO ISEC = 2, NSEC(ISURF) DY = XYZLES(2,ISEC, ISURF) - XYZLES(2,ISEC-1, ISURF) @@ -142,7 +142,13 @@ SUBROUTINE MAKESURF(ISURF) ENDDO C ELSE -C----- set spanwise spacing using overall parameters NVS(ISURF), SSPACE +C +C----- Otherwise, set spanwise spacing using the SURFACE spanwise +C parameters NVS, SSPACE +C +C This spanwise spacing is modified (fudged) to align vortex edges +C with SECTIONs as defined. This allows CONTROLs to be defined +C without bridging vortex strips C NSPACE = 2*NVS(ISURF) + 1 IF(NSPACE.GT.KPMAX) THEN @@ -177,7 +183,7 @@ SUBROUTINE MAKESURF(ISURF) IPTLOC(1) = 1 IPTLOC(NSEC(ISURF)) = NPT C -C----- fudge Glauert angles to make nodes match up exactly with interior sections +C----- fudge spacing array to make nodes match up exactly with interior sections DO ISEC = 2, NSEC(ISURF)-1 IPT1 = IPTLOC(ISEC-1) IPT2 = IPTLOC(ISEC ) @@ -187,6 +193,7 @@ SUBROUTINE MAKESURF(ISURF) STOP ENDIF C +C----- fudge spacing to this section so that nodes match up exactly with section YPT1 = YPT(IPT1) YSCALE = (YZLEN(ISEC)-YZLEN(ISEC-1)) / (YPT(IPT2)-YPT(IPT1)) DO IPT = IPT1, IPT2-1 @@ -196,6 +203,7 @@ SUBROUTINE MAKESURF(ISURF) YCP(IVS) = YZLEN(ISEC-1) + YSCALE*(YCP(IVS)-YPT1) ENDDO C +C----- check for unique spacing node for next section, if not we need more nodes IPT1 = IPTLOC(ISEC ) IPT2 = IPTLOC(ISEC+1) IF(IPT1.EQ.IPT2) THEN @@ -204,6 +212,7 @@ SUBROUTINE MAKESURF(ISURF) STOP ENDIF C +C----- fudge spacing to this section so that nodes match up exactly with section YPT1 = YPT(IPT1) YSCALE = (YPT(IPT2)-YZLEN(ISEC)) / (YPT(IPT2)-YPT(IPT1)) DO IPT = IPT1, IPT2-1 @@ -220,6 +229,19 @@ SUBROUTINE MAKESURF(ISURF) ENDDO C ENDIF +cc#ifdef USE_CPOML +C... store section counters + IF (ISURF .EQ. 1) THEN + ICNTFRST(ISURF) = 1 + ELSE + ICNTFRST(ISURF) = ICNTFRST(ISURF-1) + NCNTSEC(ISURF-1) + ENDIF + NCNTSEC(ISURF) = NSEC(ISURF) + DO ISEC = 1, NSEC(ISURF) + II = ICNTFRST(ISURF) + (ISEC-1) + ICNTSEC(II) = IPTLOC(ISEC) + ENDDO +cc#endif C C C==================================================== @@ -349,6 +371,15 @@ SUBROUTINE MAKESURF(ISURF) TANLE(idx_strip) = (XYZLER(1)-XYZLEL(1))/WIDTH TANTE(idx_strip) = (XYZLER(1)+CHORDR - XYZLEL(1)-CHORDL)/WIDTH C +cc#ifdef USE_CPOML + CHSIN = CHSINL + F1*(CHSINR-CHSINL) + CHCOS = CHCOSL + F1*(CHCOSR-CHCOSL) + AINC1(idx_strip) = ATAN2(CHSIN,CHCOS) + CHSIN = CHSINL + F2*(CHSINR-CHSINL) + CHCOS = CHCOSL + F2*(CHCOSR-CHCOSL) + AINC2(idx_strip) = ATAN2(CHSIN,CHCOS) +C +cc#endif CHSIN = CHSINL + FC*(CHSINR-CHSINL) CHCOS = CHCOSL + FC*(CHCOSR-CHCOSL) AINC(idx_strip) = ATAN2(CHSIN,CHCOS) @@ -475,6 +506,8 @@ SUBROUTINE MAKESURF(ISURF) idx_vor = IJFRST(idx_strip) DO 1505 IVC = 1, NVC(ISURF) ! NVOR = NVOR + 1 + ! change all NVOR indices into idx_vor + ! change all NSTRIP indices into idx_strip C RV1(1,idx_vor) = RLE1(1,idx_strip) & + XVR(IVC)*CHORD1(idx_strip) @@ -535,8 +568,35 @@ SUBROUTINE MAKESURF(ISURF) C C---------- TE control point used only if surface sheds a wake LVNC(idx_vor) = LFWAKE(ISURF) - idx_vor = idx_vor + 1 +C +cc#ifdef USE_CPOML +C... nodal grid associated with vortex strip (aft-panel nodes) +C... NOTE: airfoil in plane of wing, but not rotated perpendicular to dihedral; +C... retained in (x,z) plane at this point + CALL AKIMA( XLASEC(1,ISEC,ISURF), ZLASEC(1,ISEC,ISURF), NSL, + & XPT(IVC+1), ZL, DSDX ) + CALL AKIMA( XUASEC(1,ISEC,ISURF), ZUASEC(1,ISEC,ISURF), NSL, + & XPT(IVC+1), ZU, DSDX ) +C + XYN1(1,idx_vor) = RLE1(1,idx_strip) + + & XPT(IVC+1)*CHORD1(idx_strip) + XYN1(2,idx_vor) = RLE1(2,idx_strip) + ZLON1(idx_vor) = RLE1(3,idx_strip) + ZL*CHORD1(idx_strip) + ZUPN1(idx_vor) = RLE1(3,idx_strip) + ZU*CHORD1(idx_strip) +C + CALL AKIMA( XLASEC(1,ISEC+1,ISURF), ZLASEC(1,ISEC+1,ISURF), + & NSL, XPT(IVC+1), ZL, DSDX ) + CALL AKIMA( XUASEC(1,ISEC+1,ISURF), ZUASEC(1,ISEC+1,ISURF), + & NSL, XPT(IVC+1), ZU, DSDX ) + XYN2(1,idx_vor) = RLE2(1,idx_strip) + + & XPT(IVC+1)*CHORD2(idx_strip) + XYN2(2,idx_vor) = RLE2(2,idx_strip) + ZLON2(idx_vor) = RLE2(3,idx_strip) + ZL*CHORD2(idx_strip) + ZUPN2(idx_vor) = RLE2(3,idx_strip) + ZU*CHORD2(idx_strip) +C +cc#endif + idx_vor = idx_vor + 1 1505 CONTINUE C idx_strip = idx_strip + 1 @@ -723,7 +783,8 @@ SUBROUTINE SDUPL(NN, Ypt,MSG) C NNI = NN + 1 IF(NNI.GT.NFMAX) THEN - WRITE(*,*) 'SDUPL: Surface array overflow. Increase NFMAX.' + WRITE(*,*) 'SDUPL: Surface array overflow. Increase NFMAX', + & ' currently ',NFMAX STOP ENDIF C @@ -744,6 +805,7 @@ SUBROUTINE SDUPL(NN, Ypt,MSG) LFWAKE(NNI) = LFWAKE(NN) LFALBE(NNI) = LFALBE(NN) LFLOAD(NNI) = LFLOAD(NN) + LRANGE(NNI) = LRANGE(NN) LSURFSPACING(NNI) = LSURFSPACING(NN) C---- accumulate stuff for new image surface @@ -765,6 +827,16 @@ SUBROUTINE SDUPL(NN, Ypt,MSG) C C--- Image flag reversed (set to -IMAGS) for imaged surfaces IMAGS(NNI) = -IMAGS(NN) +C +cc#ifdef USE_CPOML + ICNTFRST(NNI) = ICNTFRST(NN) + NCNTSEC(NN) + NCNTSEC(NNI) = NCNTSEC(NN) + DO ISEC = 1, NCNTSEC(NNI) + IDUP = ICNTFRST(NNI) + (ISEC-1) + IORG = ICNTFRST(NN ) + (ISEC-1) + ICNTSEC(IDUP) = ICNTSEC(IORG) + ENDDO +cc#endif C YOFF = 2.0*Ypt C @@ -775,7 +847,8 @@ SUBROUTINE SDUPL(NN, Ypt,MSG) DO 100 IVS = 1, NVS(NNI) ! NSTRIP = NSTRIP + 1 IF(idx_strip.GT.NSMAX) THEN - WRITE(*,*) 'SDUPL: Strip array overflow. Increase NSMAX.' + WRITE(*,*) 'SDUPL: Strip array overflow. Increase NSMAX', + & ' currently ',NSMAX STOP ENDIF C @@ -797,6 +870,11 @@ SUBROUTINE SDUPL(NN, Ypt,MSG) TANLE(JJI) = -TANLE(JJ) AINC (JJI) = AINC(JJ) C +cc#ifdef USE_CPOML + AINC1(JJI) = AINC2(JJ) + AINC2(JJI) = AINC1(JJ) +C +cc#endif NSURFS(idx_strip) = NNI C DO N = 1, NDESIGN @@ -831,7 +909,8 @@ SUBROUTINE SDUPL(NN, Ypt,MSG) DO 80 IVC = 1, NVC(NNI) ! NVOR = NVOR + 1 IF(idx_vor.GT.NVMAX) THEN - WRITE(*,*) 'SDUPL: Vortex array overflow. Increase NVMAX.' + WRITE(*,*) 'SDUPL: Vortex array overflow. Increase NVMAX', + & ' currently ',NVMAX STOP ENDIF C @@ -862,10 +941,23 @@ SUBROUTINE SDUPL(NN, Ypt,MSG) RSGN = VREFL(JJ,N) DCONTROL(III,N) = -DCONTROL(II,N)*RSGN ENDDO +C +cc#ifdef USE_CPOML +C... nodal grid associated with vortex strip + XYN1(1,III) = XYN2(1,II) + XYN1(2,III) = -XYN2(2,II) + YOFF + XYN2(1,III) = XYN1(1,II) + XYN2(2,III) = -XYN1(2,II) + YOFF +C + ZLON1(III) = ZLON2(II) + ZUPN1(III) = ZUPN2(II) + ZLON2(III) = ZLON1(II) + ZUPN2(III) = ZUPN1(II) +cc#endif idx_vor = idx_vor + 1 - C 80 CONTINUE + idx_strip = idx_strip + 1 C 100 CONTINUE C @@ -890,7 +982,8 @@ SUBROUTINE BDUPL(NN,Ypt,MSG) C NNI = NBODY + 1 IF(NNI.GT.NBMAX) THEN - WRITE(*,*) 'BDUPL: Body array overflow. Increase NBMAX.' + WRITE(*,*) 'BDUPL: Body array overflow. Increase NBMAX', + & ' currently ',NBMAX STOP ENDIF C diff --git a/src/amode.f b/src/amode.f index 925a695..0fd9b51 100644 --- a/src/amode.f +++ b/src/amode.f @@ -29,8 +29,8 @@ SUBROUTINE MODE CHARACTER*1 ITEM, ANS, CHKEY CHARACTER*2 OPT, CHPLOT CHARACTER*4 COMAND, ITEMC - CHARACTER*120 FNOUT, FNNEW, FNSYS - CHARACTER*120 LINE, FNVB, COMARG, PROMPT, RTNEW + CHARACTER*256 FNOUT, FNNEW, FNSYS, FNVB + CHARACTER*120 LINE, COMARG, PROMPT, RTNEW CHARACTER SATYPE*50, ROTTYPE*50 C LOGICAL LPROOT(JEMAX,NRMAX), @@ -81,6 +81,11 @@ SUBROUTINE MODE c ROBINV = 0. C OVERLAY = .FALSE. + + if (LVERBOSE) then + write(*,*) 3, parnam(ipixx), ' ', parunch(ipixx) + endif + C C================================================================= C---- start of user interaction loop @@ -95,7 +100,7 @@ SUBROUTINE MODE C C C - WRITE(*,1052) + WRITE(*,1052) LSVMOV 1052 FORMAT( & ' ==========================================================' & //' "#" select run case for eigenmode analysis (0 = all)' @@ -108,6 +113,7 @@ SUBROUTINE MODE & //' A nnotate current plot' & /' H ardcopy current plot' & /' T ime-integration parameters' + & /' G enerate hardcopy movie toggle', L3 & //' S ystem matrix output' & /' W rite eigenvalues to file' & /' D ata file overlay toggle' @@ -115,7 +121,7 @@ SUBROUTINE MODE & /' U nzoom') C C A B C D E F G H I J K L M N O P Q R S T U V W X Y Z -C x x x x x x x x x x x x x x +C x x x x x x x x x x x x x x x x 810 CONTINUE CALL ASKC(' .MODE^',COMAND,COMARG) @@ -205,7 +211,7 @@ SUBROUTINE MODE ENDIF C NITER = 10 - INFO = 0 + INFO = 1 CALL EXEC(NITER,INFO,IR) C IF(COMAND .EQ. 'N ') THEN @@ -463,6 +469,10 @@ SUBROUTINE MODE ENDIF GO TO 70 C +C----------------------------------------------------------------------- + ELSEIF(COMAND.EQ.'G ') THEN + LSVMOV = .NOT. LSVMOV +C C------------------------------------------------------------------- C---- write system matrices ELSEIF(COMAND.EQ.'S ' .OR. @@ -489,7 +499,7 @@ SUBROUTINE MODE CDREF = PARVAL(IPCD0,IR) C NITER = 10 - INFO = 0 + INFO = 1 CALL EXEC(NITER,INFO,IR) C IF(COMAND.EQ.'S ') THEN @@ -528,25 +538,25 @@ SUBROUTINE MODE C C------------------------------------------------------------------- C---- write eigenvalues -C ELSEIF(COMAND.EQ.'W ') THEN -C IF(FEVDEF(1:1).EQ.' ') THEN -C C------ set default filename -C KDOT = INDEX(FILDEF,'.') -C IF(KDOT.EQ.0) THEN -C CALL SLEN(FILDEF,NFIL) -C FEVDEF = FILDEF(1:NFIL) // '.eig' -C ELSE -C FEVDEF = FILDEF(1:KDOT) // 'eig' -C ENDIF -C ENDIF -C C -C CALL SLEN(FEVDEF,NFE) -C NFE = MAX( NFE , 1 ) -C WRITE(*,2040) FEVDEF(1:NFE) -C 2040 FORMAT(' Enter eigenvalue save filename: ', A) -C READ (*,1000) FNNEW -C IF(FNNEW.NE.' ') FEVDEF = FNNEW -C CALL EIGOUT(FEVDEF, IRUN1,IRUN2, SAVED) + ELSEIF(COMAND.EQ.'W ') THEN + IF(FEVDEF(1:1).EQ.' ') THEN +C------ set default filename + KDOT = INDEX(FILDEF,'.') + IF(KDOT.EQ.0) THEN + CALL SLEN(FILDEF,NFIL) + FEVDEF = FILDEF(1:NFIL) // '.eig' + ELSE + FEVDEF = FILDEF(1:KDOT) // 'eig' + ENDIF + ENDIF +C + CALL SLEN(FEVDEF,NFE) + NFE = MAX( NFE , 1 ) + WRITE(*,2040) FEVDEF(1:NFE) + 2040 FORMAT(' Enter eigenvalue save filename: ', A) + READ (*,1000) FNNEW + IF(FNNEW.NE.' ') FEVDEF = FNNEW + CALL EIGOUT(FEVDEF, IRUN1,IRUN2, SAVED) C C------------------------------------------------------------------- C C---- toggle overlay flag @@ -988,7 +998,7 @@ SUBROUTINE SYSMAT(IR,ASYS,BSYS,RSYS,NSYS) & + RIINV(K,2)*WXH_U(2,IU) & + RIINV(K,3)*WXH_U(3,IU) ENDDO -C + DO N = 1, NCONTROL MIF_D(K,N) = & MAINV(K,1)*CXTOT_D(N)*QS @@ -1043,6 +1053,11 @@ SUBROUTINE SYSMAT(IR,ASYS,BSYS,RSYS,NSYS) DO N = 1, NCONTROL BSYS(IEQ,N) = MIF_D(K,N) ENDDO + +c MIF_U(1,2) +c write(*,*) +c write(*,*) '* 1u', MIF_U(K,2), PRF_U(K,2) + C C---- y-acceleration IEQ = JEV @@ -1234,7 +1249,6 @@ SUBROUTINE SYSMAT(IR,ASYS,BSYS,RSYS,NSYS) & + TT_ANG(K,2,3)*VINF(2) & + TT_ANG(K,3,3)*VINF(3) )*VEE C - c write(*,*) 'H ', H c write(*,*) 'WxH', WXH c write(*,*) 'I-1 M ', RIM @@ -1250,7 +1264,9 @@ SUBROUTINE SYSMAT(IR,ASYS,BSYS,RSYS,NSYS) SUBROUTINE APPMAT(IR,ASYS,BSYS,RSYS,NSYS) C------------------------------------------------------------------ -C Computes system matrices for run case IR. +C Computes system matrices for run case IR, +C using the approximate expression in Etkin. +C C Current forces and derivatives are assumed to be correct. C------------------------------------------------------------------ INCLUDE 'AVL.INC' @@ -1526,11 +1542,22 @@ SUBROUTINE APPMAT(IR,ASYS,BSYS,RSYS,NSYS) SUBROUTINE SYSSHO(LU,ASYS,BSYS,RSYS,NSYS) C------------------------------------------------------------------ -C Computes eigenvalues and eigenvectors for run case IR. -C Current forces and derivatives are assumed to be correct. +C Prints out state-system matrices "A" and "B" +C In an organized manner. C------------------------------------------------------------------ INCLUDE 'AVL.INC' REAL*8 ASYS(JEMAX,JEMAX),BSYS(JEMAX,NDMAX),RSYS(JEMAX) + REAL*8 USGN(JEMAX) +C + DO I = 1, NSYS + USGN(I) = 1.0 + ENDDO + USGN(JEU) = -1.0 + USGN(JEW) = -1.0 + USGN(JEP) = -1.0 + USGN(JER) = -1.0 + USGN(JEX) = -1.0 + USGN(JEZ) = -1.0 C WRITE(LU,*) WRITE(LU,1100) @@ -1542,8 +1569,10 @@ SUBROUTINE SYSSHO(LU,ASYS,BSYS,RSYS,NSYS) 1100 FORMAT(1X,A,A,A,1X,'|',2X,12A12) C DO I = 1, NSYS - WRITE(LU,1200) (ASYS(I,J), J=1, NSYS), - & (BSYS(I,N), N=1, NCONTROL) + WRITE(LU,1200) + & (ASYS(I,J)*USGN(I)*USGN(J), J=1, NSYS), + & (BSYS(I,N)*USGN(I) , N=1, NCONTROL) +c & , RSYS(I)*USGN(I) 1200 FORMAT(1X,12F10.4,3X,12G12.4) ENDDO C diff --git a/src/aoml.f b/src/aoml.f new file mode 100644 index 0000000..b89a46b --- /dev/null +++ b/src/aoml.f @@ -0,0 +1,470 @@ +C*********************************************************************** +C Module: aoml.f +C +C Copyright (C) 2002 Mark Drela, Harold Youngren +C +C This program is free software; you can redistribute it and/or modify +C it under the terms of the GNU General Public License as published by +C the Free Software Foundation; either version 2 of the License, or +C (at your option) any later version. +C +C This program is distributed in the hope that it will be useful, +C but WITHOUT ANY WARRANTY; without even the implied warranty of +C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +C GNU General Public License for more details. +C +C You should have received a copy of the GNU General Public License +C along with this program; if not, write to the Free Software +C Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. +C*********************************************************************** + + + + SUBROUTINE CPOML +C + INCLUDE 'AVL.INC' +C + LOGICAL LRANGEALL +C +C... check that all surfaces use full range of airfoil definitions + LRANGEALL = .TRUE. + DO ISURF = 1, NSURF + LRANGEALL = LRANGEALL .AND. LRANGE(ISURF) + ENDDO + IF (.NOT.LRANGEALL) THEN + WRITE(*,*) 'ERROR in CPOML: implemented only for surfaces ', + & 'defined using full range of input airfoils' + WRITE(*,*) ' returning without writing OML output' + RETURN + ENDIF +C + CALL CPTHK() + CALL CPDUMP() +C + RETURN + END + + + SUBROUTINE CPTHK +C +C...PURPOSE compute sectional pressure coefficients due to thickness +C via constant-source panels along thickness surface +C NOTE: airfoil assumed in (x,z) plane even if there is dihedral +C +C...INPUT +C +C...OUTPUT CPT thickness based pressure coefficient (in common) +C +C...COMMENTS +C + INCLUDE 'AVL.INC' +C + PARAMETER (NCMAX=256) + DIMENSION AICT(NCMAX,NCMAX) ! AIC for no normal flow BC + DIMENSION BICT(NCMAX,NCMAX) ! AIC for tangential flow (used after solve for Cp) + DIMENSION RHS_T(NCMAX), WORK(NCMAX) + DIMENSION QSINF(NCMAX) ! freestream tangential velocity + DIMENSION SRCTHK(NCMAX) +C + DO ISTRIP = 1, NSTRIP +C + I1 = IJFRST(ISTRIP) + NVC_strp = NVSTRP(ISTRIP) + IF (NVC_strp.GT.NCMAX) THEN + WRITE(*,*) '* CPTHK: Array overflow. Increase NCMAX to', NVC_strp + STOP + ENDIF +C + XLE = 0.5*(RLE1(1,ISTRIP) + RLE2(1,ISTRIP)) + ZLE = 0.5*(RLE1(3,ISTRIP) + RLE2(3,ISTRIP)) +C + DO II = 1, NVC_strp + I = I1 + (II-1) +C +C... control point: panel center + IF (II.EQ.1) THEN + X = 0.5*(XLE + 0.5*(XYN1(1,I) + XYN2(1,I))) - XLE + ZLO = 0.5*(ZLE + 0.5*(ZLON1(I) + ZLON2(I))) - ZLE + ZUP = 0.5*(ZLE + 0.5*(ZUPN1(I) + ZUPN2(I))) - ZLE + ELSE + X = 0.25*(XYN1(1,I-1) + XYN2(1,I-1) + XYN1(1,I) + XYN2(1,I)) + ZLO = 0.25*(ZLON1(I-1) + ZLON2(I-1) + ZLON1(I) + ZLON2(I)) + ZUP = 0.25*(ZUPN1(I-1) + ZUPN2(I-1) + ZUPN1(I) + ZUPN2(I)) + X = X - XLE + ZLO = ZLO - ZLE + ZUP = ZUP - ZLE + ENDIF + Z = 0.5*(ZUP - ZLO) +C +C... unit-normal to thickness at panel center + IF (II.EQ.1) THEN + DX = 0.5*(XYN1(1,I) + XYN2(1,I)) - XLE + DZLO = 0.5*(ZLON1(I) + ZLON2(I)) - ZLE + DZUP = 0.5*(ZUPN1(I) + ZUPN2(I)) - ZLE + ELSE + DX = 0.5*(XYN1(1,I) - XYN1(1,I-1) + XYN2(1,I) - XYN2(1,I-1)) + DZLO = 0.5*(ZLON1(I) - ZLON1(I-1) + ZLON2(I) - ZLON2(I-1)) + DZUP = 0.5*(ZUPN1(I) - ZUPN1(I-1) + ZUPN2(I) - ZUPN2(I-1)) + ENDIF + DZ = 0.5*(DZUP - DZLO) + + DEN = SQRT(DX*DX + DZ*DZ) + ESX = DX/DEN + ESZ = DZ/DEN + ENX = -ESZ + ENZ = ESX +C +C... thickness solution at zero alpha + VINF1 = 1 + VINF3 = 0 + RHS_T(II) = ENX*VINF1 + ENZ*VINF3 + QSINF(II) = ESX*VINF1 + ESZ*VINF3 ! needed for Cp +C + DO JJ = 1, NVC_strp + J = I1 + (JJ-1) +C + IF (JJ.EQ.1) THEN + X1 = 0 + Z1 = 0 + ELSE + X1 = 0.5*(XYN1(1,J-1) + XYN2(1,J-1)) - XLE + Z1LO = 0.5*(ZLON1(J-1) + ZLON2(J-1)) - ZLE + Z1UP = 0.5*(ZUPN1(J-1) + ZUPN2(J-1)) - ZLE + Z1 = 0.5*(Z1UP - Z1LO) + ENDIF +C + X2 = 0.5*(XYN1(1,J) + XYN2(1,J)) - XLE + Z2LO = 0.5*(ZLON1(J) + ZLON2(J)) - ZLE + Z2UP = 0.5*(ZUPN1(J) + ZUPN2(J)) - ZLE + Z2 = 0.5*(Z2UP - Z2LO) + IF (JJ.EQ.II) THEN + U = 0.5*ENX; !! trick: needed to ensure proper velocity at ii=jj panel; + W = 0.5*ENZ; !! if (x,z) below panel, then get -0.5 (which is bad) + ELSE + CALL SRCPANEL( X, Z, X1, Z1, X2, Z2, U, W ) + ENDIF + AICT(II,JJ) = ENX*U + ENZ*W + BICT(II,JJ) = ESX*U + ESZ*W +C + CALL SRCPANEL( X, Z, X1, -Z1, X2, -Z2, U, W ) + AICT(II,JJ) = AICT(II,JJ) + ENX*U + ENZ*W + BICT(II,JJ) = BICT(II,JJ) + ESX*U + ESZ*W + ENDDO + ENDDO +C + CALL LUDCMP(NCMAX,NVC_strp,AICT,IAPIV,WORK) + CALL BAKSUB(NCMAX,NVC_strp,AICT,IAPIV,RHS_T) +C + DO II = 1, NVC_strp + SRCTHK(II) = -RHS_T(II) + ENDDO +C + DO II = 1, NVC_strp + I = I1 + (II-1) +C + QS = QSINF(II) + DO JJ = 1, NVC_strp + QS = QS + BICT(II,JJ)*SRCTHK(JJ) + ENDDO + CPT(I) = 1 - QS*QS + ENDDO + ENDDO ! ISTRIP +C + RETURN + END ! CPTHK + + + SUBROUTINE CPDUMP +C +C...PURPOSE output OML upper and lower grid and pressure coefficients +C +C...INPUT +C +C...OUTPUT +C +C...COMMENTS C code reader in read_cpoml.c +C + INCLUDE 'AVL.INC' +C + LU = 12 + OPEN(LU, FILE='cpoml.dat', FORM='FORMATTED', STATUS='UNKNOWN') +C + WRITE(LU,'(A)') 'CPOML' +C + WRITE(LU,'(A)') 'VERSION 1.0' + WRITE(LU,'(I6,6X,A)') NSURF, ' | surfaces' + DO ISURF = 1, NSURF + NVC_chord = NK(ISURF) +C +C... determine if surface is L-to-R or R-to-L + IF (NJ(ISURF) .EQ. 1) THEN + ISTRIP0 = JFRST(ISURF) + ISTRIP1 = ISTRIP0 + NJ(ISURF) - 1 + ISTEP = 1 + ELSE + ISTRIP0 = JFRST(ISURF) + ISTRIP1 = ISTRIP0 + NJ(ISURF) - 1 +C + Y0 = RLE1(2,ISTRIP0) + Y1 = RLE1(2,ISTRIP1) + IF (Y1 .GT. Y0) THEN + ISTEP = 1 ! L-to-R + ELSE + ISTRIP1 = JFRST(ISURF) + ISTRIP0 = ISTRIP1 + NJ(ISURF) - 1 + ISTEP = -1 ! R-to-L + ENDIF + ENDIF +C + WRITE(LU,'(A)') 'SURFACE' + WRITE(LU,'(A)') STITLE(ISURF) + WRITE(LU,'(I6,6X,A)') LSCOMP(ISURF), ' | component' + WRITE(LU,'(2I6,A)') NJ(ISURF), NK(ISURF), + & ' | elements (span, chord)' + WRITE(LU,'(I6,6X,A)') IMAGS(ISURF), + & ' | imags (<0 if Y-duplicated)' +C + NSEC_surf = NCNTSEC(ISURF) + II = ICNTFRST(ISURF) + WRITE(LU,'(I6,6X,A)') NSEC_surf, ' | section indices' + IF (ISTEP .GT. 0) THEN + WRITE(LU,*) (ICNTSEC(II + ISEC-1), ISEC=1,NSEC_surf) + ELSE + NSPAN = NJ(ISURF) + 1 + WRITE(LU,*) (NSPAN - ICNTSEC(II + ISEC-1) + 1, + & ISEC=NSEC_surf,1,-1) + ENDIF +C + WRITE(LU,'(A,1X,A)') 'VERTEX_GRID', + & '(x_lo, x_up, y_lo, y_up, z_lo, z_up)' + DO J = ISTRIP0, ISTRIP1, ISTEP + I1 = IJFRST(J) +C + DYLE = RLE2(2,J) - RLE1(2,J) + DZLE = RLE2(3,J) - RLE1(3,J) + + CSA = COS(AINC1(J)) + SNA = SIN(AINC1(J)) +C + XLE = RLE1(1,J) + YLE = RLE1(2,J) + ZLE = RLE1(3,J) + WRITE(LU,'(6(ES23.15))') XLE, XLE, YLE, YLE, ZLE, ZLE +C + DO II = 1, NVC_chord + I = I1 + (II-1) +C + X0 = XYN1(1,I) + Y0 = XYN1(2,I) + ZLO0 = ZLON1(I) + ZUP0 = ZUPN1(I) +C +C... rotate airfoil in (y,z) so that it is perpendicular to dihedral + DY = XYN2(2,I) - XYN1(2,I) + DZ = 0.5*(ZLON2(I) - ZLON1(I)) + & + 0.5*(ZUPN2(I) - ZUPN1(I)) + CSD = DY/SQRT(DY*DY + DZ*DZ) + SND = DZ/SQRT(DY*DY + DZ*DZ) + YLOD = YLE + (Y0 - YLE)*CSD - (ZLO0 - ZLE)*SND + YUPD = YLE + (Y0 - YLE)*CSD - (ZUP0 - ZLE)*SND + ZLOD = ZLE - (Y0 - YLE)*SND + (ZLO0 - ZLE)*CSD + ZUPD = ZLE - (Y0 - YLE)*SND + (ZUP0 - ZLE)*CSD +C +C... rotate airfoil in (x,z) for twist + XLO = XLE + (X0 - XLE)*CSA + (ZLOD - ZLE)*SNA + XUP = XLE + (X0 - XLE)*CSA + (ZUPD - ZLE)*SNA + ZLO = ZLE - (X0 - XLE)*SNA + (ZLOD - ZLE)*CSA + ZUP = ZLE - (X0 - XLE)*SNA + (ZUPD - ZLE)*CSA + YLO = YLOD + YUP = YUPD + WRITE(LU,'(6(ES23.15))') XLO, XUP, YLO, YUP, ZLO, ZUP + ENDDO + ENDDO +C + J = ISTRIP1 + I1 = IJFRST(J) +C + CSA = COS(AINC2(J)) + SNA = SIN(AINC2(J)) +C + XLE = RLE2(1,J) + YLE = RLE2(2,J) + ZLE = RLE2(3,J) + WRITE(LU,'(6(ES23.15))') XLE, XLE, YLE, YLE, ZLE, ZLE +C + DO II = 1, NVC_chord + I = I1 + (II-1) +C + X0 = XYN2(1,I) + Y0 = XYN2(2,I) + ZLO0 = ZLON2(I) + ZUP0 = ZUPN2(I) +C +C... rotate airfoil in (y,z) so that it is perpendicular to dihedral + DY = XYN2(2,I) - XYN1(2,I) + DZ = 0.5*(ZLON2(I) - ZLON1(I)) + & + 0.5*(ZUPN2(I) - ZUPN1(I)) + CSD = DY/SQRT(DY*DY + DZ*DZ) + SND = DZ/SQRT(DY*DY + DZ*DZ) + YLOD = YLE + (Y0 - YLE)*CSD - (ZLO0 - ZLE)*SND + YUPD = YLE + (Y0 - YLE)*CSD - (ZUP0 - ZLE)*SND + ZLOD = ZLE - (Y0 - YLE)*SND + (ZLO0 - ZLE)*CSD + ZUPD = ZLE - (Y0 - YLE)*SND + (ZUP0 - ZLE)*CSD +C +C... rotate airfoil in (x,z) for twist + XLO = XLE + (X0 - XLE)*CSA + (ZLOD - ZLE)*SNA + XUP = XLE + (X0 - XLE)*CSA + (ZUPD - ZLE)*SNA + ZLO = ZLE - (X0 - XLE)*SNA + (ZLOD - ZLE)*CSA + ZUP = ZLE - (X0 - XLE)*SNA + (ZUPD - ZLE)*CSA + YLO = YLOD + YUP = YUPD + WRITE(LU,'(6(ES23.15))') XLO, XUP, YLO, YUP, ZLO, ZUP + ENDDO +C + WRITE(LU,'(A,1X,A)') 'ELEMENT_CP', + & '(x_lo, x_up, y_lo, y_up, z_lo, z_up, cp_lo, cp_up)' + DO J = ISTRIP0, ISTRIP1, ISTEP + I1 = IJFRST(J) +C + CSA = COS(AINC(J)) + SNA = SIN(AINC(J)) +C + XLE = 0.5*(RLE1(1,J) + RLE2(1,J)) + YLE = 0.5*(RLE1(2,J) + RLE2(2,J)) + ZLE = 0.5*(RLE1(3,J) + RLE2(3,J)) +C + X0 = 0.5*(XLE + 0.5*(XYN1(1,I1) + XYN2(1,I1))) + Y0 = 0.5*(YLE + 0.5*(XYN1(2,I1) + XYN2(2,I1))) + ZLO0 = 0.5*(ZLE + 0.5*(ZLON1(I1) + ZLON2(I1))) + ZUP0 = 0.5*(ZLE + 0.5*(ZUPN1(I1) + ZUPN2(I1))) +C +C... rotate airfoil in (y,z) so that it is perpendicular to dihedral + DY = XYN2(2,I1) - XYN1(2,I1) + DZ = 0.5*(ZLON2(I1) - ZLON1(I1)) + & + 0.5*(ZUPN2(I1) - ZUPN1(I1)) + CSD = DY/SQRT(DY*DY + DZ*DZ) + SND = DZ/SQRT(DY*DY + DZ*DZ) + YLOD = YLE + (Y0 - YLE)*CSD - (ZLO0 - ZLE)*SND + YUPD = YLE + (Y0 - YLE)*CSD - (ZUP0 - ZLE)*SND + ZLOD = ZLE - (Y0 - YLE)*SND + (ZLO0 - ZLE)*CSD + ZUPD = ZLE - (Y0 - YLE)*SND + (ZUP0 - ZLE)*CSD +C +C... rotate airfoil in (x,z) for twist + XLO = XLE + (X0 - XLE)*CSA + (ZLO0 - ZLE)*SNA + XUP = XLE + (X0 - XLE)*CSA + (ZUP0 - ZLE)*SNA + ZLO = ZLE - (X0 - XLE)*SNA + (ZLO0 - ZLE)*CSA + ZUP = ZLE - (X0 - XLE)*SNA + (ZUP0 - ZLE)*CSA + YLO = YLOD + YUP = YUPD +C + CPU = CPT(I1) + 0.5*DCP(I1) + CPL = CPT(I1) - 0.5*DCP(I1) + WRITE(LU,'(8(ES23.15))') XLO,XUP, YLO,YUP, ZLO,ZUP, CPL,CPU +C + DO II = 2, NVC_chord + I = I1 + (II-1) +C + X0 = 0.5*(XYN1(1,I-1) + XYN2(1,I-1)) + Y0 = 0.5*(XYN1(2,I-1) + XYN2(2,I-1)) + X1 = 0.5*(XYN1(1,I ) + XYN2(1,I )) + Y1 = 0.5*(XYN1(2,I ) + XYN2(2,I )) + X0 = 0.5*(X0 + X1) + Y0 = 0.5*(Y0 + Y1) +C + ZLO0 = 0.5*(ZLON1(I-1) + ZLON2(I-1)) + ZUP0 = 0.5*(ZUPN1(I-1) + ZUPN2(I-1)) + ZLO1 = 0.5*(ZLON1(I ) + ZLON2(I )) + ZUP1 = 0.5*(ZUPN1(I ) + ZUPN2(I )) + ZLO0 = 0.5*(ZLO0 + ZLO1) + ZUP0 = 0.5*(ZUP0 + ZUP1) +C +C... rotate airfoil in (y,z) so that it is perpendicular to dihedral + DY = XYN2(2,I) - XYN1(2,I) + DZ = 0.5*(ZLON2(I) - ZLON1(I)) + & + 0.5*(ZUPN2(I) - ZUPN1(I)) + CSD = DY/SQRT(DY*DY + DZ*DZ) + SND = DZ/SQRT(DY*DY + DZ*DZ) + YLOD = YLE + (Y0 - YLE)*CSD - (ZLO0 - ZLE)*SND + YUPD = YLE + (Y0 - YLE)*CSD - (ZUP0 - ZLE)*SND + ZLOD = ZLE - (Y0 - YLE)*SND + (ZLO0 - ZLE)*CSD + ZUPD = ZLE - (Y0 - YLE)*SND + (ZUP0 - ZLE)*CSD +C +C... rotate airfoil in (x,z) for twist + XLO = XLE + (X0 - XLE)*CSA + (ZLOD - ZLE)*SNA + XUP = XLE + (X0 - XLE)*CSA + (ZUPD - ZLE)*SNA + ZLO = ZLE - (X0 - XLE)*SNA + (ZLOD - ZLE)*CSA + ZUP = ZLE - (X0 - XLE)*SNA + (ZUPD - ZLE)*CSA + YLO = YLOD + YUP = YUPD +C + CPU = CPT(I) + 0.5*DCP(I) + CPL = CPT(I) - 0.5*DCP(I) + WRITE(LU,'(8(ES23.15))') XLO,XUP, YLO,YUP, ZLO,ZUP, CPL,CPU + ENDDO + ENDDO + ENDDO +C + CLOSE(LU) +C + RETURN + END ! CPDUMP + + + SUBROUTINE SRCPANEL( X, Z, X1, Z1, X2, Z2, U, W ) +C +C...PURPOSE compute velocities induced by 2D unit-strength constant-source panel +C +C...INPUT x, z control point +C x1, z1 panel endpoint 1 +C x2, z2 panel endpoint 2 +C +C...OUTPUT u, w induced velocity in global coordinates +C +C...COMMENTS offset from endpoint 1 mitigates issues with atan branch cut +C +C... rotate to local panel coordinates + DEN = SQRT((X2 - X1)**2 + (Z2 - Z1)**2) + CS = (X2 - X1)/DEN + SN = (Z2 - Z1)/DEN +C + XP = CS*(X - X1) + SN*(Z - Z1) + ZP = -SN*(X - X1) + CS*(Z - Z1) + X1P = 0 + Z1P = 0 + X2P = DEN + Z2P = 0 +C +C... velocity induced by unit-strength constant-source panel + PI = 4*ATAN(1.) + R1SQ = (XP - X1P)**2 + (ZP - Z1P)**2 + R2SQ = (XP - X2P)**2 + (ZP - Z2P)**2 + TH1 = ATAN2(ZP - Z1P, XP - X1P) + TH2 = ATAN2(ZP - Z2P, XP - X2P) + UP = LOG(R1SQ/R2SQ) / (4*PI) + WP = (TH2 - TH1) / (2*PI) +C +C... rotate velocity to global + U = CS*UP - SN*WP + W = SN*UP + CS*WP + + RETURN + END ! SRCPANEL + + + SUBROUTINE DUMPMARKET(NNDIM,NN,A) + DIMENSION A(NNDIM,NNDIM) +C + OPEN(10, FILE='jac.mtx', FORM='FORMATTED', STATUS='UNKNOWN') + WRITE(10,100) + WRITE(10,*) NN, NN, NN*NN + DO I = 1, NN + DO J = 1, NN + WRITE(10,*) I, J, A(I,J) + ENDDO + ENDDO + CLOSE(10) +C +100 FORMAT('%%MatrixMarket matrix coordinate real general') + RETURN + END + diff --git a/src/aoper.f b/src/aoper.f index 5e30231..3e7fc61 100644 --- a/src/aoper.f +++ b/src/aoper.f @@ -25,19 +25,18 @@ SUBROUTINE OPER INCLUDE 'AVL.INC' C INCLUDE 'AVLPLT.INC' LOGICAL ERROR, LCERR, LCERRI, LWRIT, LMATCH, LMATCH_TEST + LOGICAL USEMRF C - CHARACTER*1 ANS CHARACTER*2 OPT CHARACTER*4 COMAND, COMAND_TEST, ITEMC - CHARACTER*120 FNOUT, FNDER, FNNEW - CHARACTER*120 LINE, FNVB, COMARG, COMARG_TEST, CRUN, PROMPT, RTNEW + CHARACTER*256 FNOUT, FNDER, FNNEW, FNVB + CHARACTER*120 LINE, COMARG, COMARG_TEST, CRUN, PROMPT, RTNEW CHARACTER*50 SATYPE, ROTTYPE C LOGICAL LOWRIT C REAL RINPUT(20), RINP(20) - INTEGER IINPUT(20), IINP(20), IRUN_TEST C IF(.NOT.LGEO) THEN @@ -62,16 +61,18 @@ SUBROUTINE OPER C LPLOT = .FALSE. LWRIT = .FALSE. LSOL = .FALSE. +C + USEMRF = .FALSE. C FNOUT = ' ' C C================================================================= C---- start of user interaction loop 800 CONTINUE -C C +C LCERR = .FALSE. -C C -C C +C +C CALL CFRAC(IRUN,NRUN,CRUN,NPR) C if(lverbose)then @@ -79,7 +80,7 @@ SUBROUTINE OPER CALL CONLST(IRUN) end if C WRITE(*,1052) - +C 1050 FORMAT( & /' Operation of run case ',A,': ', A & /' ==========================================================') @@ -100,7 +101,8 @@ SUBROUTINE OPER & /' DE design changes FE element forces ' & /' O ptions FB body forces ' & /' HM hinge moments ' - & /' VM strip shear,moment ') + & /' VM strip shear,moment ' + & /' MRF machine-readable format CPOM OML surface pressures') C C A B C D E F G H I J K L M N O P Q R S T U V W X Y Z C x x x x x x x x x x x x x x x x x x x @@ -395,7 +397,11 @@ SUBROUTINE OPER C ELSEIF(LU.EQ.0) THEN C WRITE(*,*) '* Data not written' C ELSE -C CALL OUTTOT(LU) +C IF(USEMRF) THEN +C CALL MRFTOT(LU, 'TOT') +C ELSE +C CALL OUTTOT(LU) +C ENDIF C IF(LU.NE.5 .AND. LU.NE.6) CLOSE(LU) C ENDIF C C @@ -414,7 +420,11 @@ SUBROUTINE OPER C ELSEIF(LU.EQ.0) THEN C WRITE(*,*) '* Data not written' C ELSE +C IF(USEMRF) THEN +C CALL MRFSURF(LU) +C ELSE C CALL OUTSURF(LU) +C ENDIF C IF(LU.NE.5 .AND. LU.NE.6) CLOSE(LU) C ENDIF C C @@ -433,7 +443,11 @@ SUBROUTINE OPER C ELSEIF(LU.EQ.0) THEN C WRITE(*,*) '* Data not written' C ELSE +C IF(USEMRF) THEN +C CALL MRFSTRP(LU) +C ELSE C CALL OUTSTRP(LU) +C ENDIF C IF(LU.NE.5 .AND. LU.NE.6) CLOSE(LU) C ENDIF C C @@ -452,7 +466,11 @@ SUBROUTINE OPER C ELSEIF(LU.EQ.0) THEN C WRITE(*,*) '* Data not written' C ELSE +C IF(USEMRF) THEN +C CALL MRFELE(LU) +C ELSE C CALL OUTELE(LU) +C ENDIF C IF(LU.NE.5 .AND. LU.NE.6) CLOSE(LU) C ENDIF C C @@ -470,8 +488,12 @@ SUBROUTINE OPER C WRITE(*,*) '* Filename error *' C ELSEIF(LU.EQ.0) THEN C WRITE(*,*) '* Data not written' +C ELSE +C IF(USEMRF) THEN +C CALL MRFBODY(LU) C ELSE C CALL OUTBODY(LU) +C ENDIF C IF(LU.NE.5 .AND. LU.NE.6) CLOSE(LU) C ENDIF C C @@ -490,7 +512,11 @@ SUBROUTINE OPER C ELSEIF(LU.EQ.0) THEN C WRITE(*,*) '* Data not written' C ELSE +C IF(USEMRF) THEN +C CALL MRFHINGE(LU) +C ELSE C CALL OUTHINGE(LU) +C ENDIF C IF(LU.NE.5 .AND. LU.NE.6) CLOSE(LU) C ENDIF C C @@ -509,7 +535,11 @@ SUBROUTINE OPER C ELSEIF(LU.EQ.0) THEN C WRITE(*,*) '* Data not written' C ELSE -C CALL OUTVM(LU) +C IF(USEMRF) THEN +C CALL MRFVM(LU) +C ELSE +C CALL OUTVM(LU) +C ENDIF C IF(LU.NE.5 .AND. LU.NE.6) CLOSE(LU) C ENDIF C C @@ -528,7 +558,11 @@ SUBROUTINE OPER C ELSEIF(LU.EQ.0) THEN C WRITE(*,*) '* Data not written' C ELSE -C CALL OUTCNC(LU) +C IF(USEMRF) THEN +C CALL MRFCNC(LU) +C ELSE +C CALL OUTCNC(LU) +C ENDIF C IF(LU.NE.5 .AND. LU.NE.6) CLOSE(LU) C ENDIF C C @@ -672,13 +706,34 @@ SUBROUTINE OPER C C C ENDIF C C +C IF(USEMRF) THEN +C IF(LPTOT) CALL MRFTOT(LUOUT, 'TOT') +C IF(LPSURF) CALL MRFSURF(LUOUT) +C IF(LPSTRP) CALL MRFSTRP(LUOUT) +C IF(LPELE) CALL MRFELE(LUOUT) +C ELSE C IF(LPTOT) CALL OUTTOT(LUOUT) C IF(LPSURF) CALL OUTSURF(LUOUT) C IF(LPSTRP) CALL OUTSTRP(LUOUT) C IF(LPELE) CALL OUTELE(LUOUT) C ccc IF(LPDERIV) CALL DERMAT(LUOUT) -C Cs +C ENDIF +C C C C------------------------------------------------------ +C ELSE IF(COMAND.EQ.'MRF ') THEN +C C------ all subsequent data file writes are full-precision machine readable format +C USEMRF = .TRUE. +C WRITE(*,'(A,A)') +C & 'MRF: all subsequent output in full-precision ', +C & 'machine readable format' +C cc#ifdef USE_CPOML +C C +C C------------------------------------------------------ +C ELSE IF(COMAND.EQ.'CPOM') THEN +C C------ write OML surface pressures to a file +C CALL CPOML() +C cc#endif +CC C------------------------------------------------------ C ELSE IF(COMAND.EQ.'RE ') THEN C C------ Change reference data C 89 WRITE(*,2090) SREF,CREF,BREF @@ -745,7 +800,7 @@ SUBROUTINE OPER C GO TO 800 C RETURN - END + END ! OPER @@ -778,7 +833,7 @@ SUBROUTINE CONLST(IR) & ' ',A,' -> ', A, '=', G12.4, 1X, A) 1030 FORMAT( & ' ------------ ------------------------') - END + END ! CONLST @@ -902,7 +957,7 @@ SUBROUTINE CONSET(COMAND,COMARG) !,LMATCH,IR) - SUBROUTINE EXEC(NITER,INFO,IR) + SUBROUTINE EXEC(NITER,INFO,IR) C--------------------------------------------------- C Solves for the flow condition specified by C the global operating parameters: @@ -924,7 +979,9 @@ SUBROUTINE EXEC(NITER,INFO,IR) C C---- convergence epsilon, max angle limit (radians) - DATA EPS, DMAX / 0.00002, 1.0 / + ! EPS replaced by EXEC_TOL + ! DATA EPS, DMAX / 0.00002, 1.0 / + DATA DMAX / 1.0 / C IF(LNASA_SA) THEN @@ -1091,6 +1148,7 @@ SUBROUTINE EXEC(NITER,INFO,IR) VSYS(IV,IVROTX) = CLTOT_U(4) VSYS(IV,IVROTY) = CLTOT_U(5) VSYS(IV,IVROTZ) = CLTOT_U(6) + C DO N = 1, NCONTROL NV = IVTOT + N @@ -1137,6 +1195,7 @@ SUBROUTINE EXEC(NITER,INFO,IR) C DO N = 1, NCONTROL NV = IVTOT + N + ! no *DIR here because we do *DIR in AERO now VSYS(IV,NV) = (CRTOT_D(N)*CA + CNTOT_D(N)*SA) ENDDO C @@ -1180,6 +1239,7 @@ SUBROUTINE EXEC(NITER,INFO,IR) C DO N = 1, NCONTROL NV = IVTOT + N + ! no *DIR here because we do *DIR in AERO now VSYS(IV,NV) = (CNTOT_D(N)*CA - CRTOT_D(N)*SA) ENDDO C @@ -1200,6 +1260,12 @@ SUBROUTINE EXEC(NITER,INFO,IR) C 100 CONTINUE C + +c write(*,*) +c do k = 1, nvtot +c write(*,'(1x,40f9.4)') (vsys(k,l), l=1, nvtot), vres(k) +c enddo +c write(*,*) C C------ LU-factor, and back-substitute RHS CALL LUDCMP(IVMAX,NVTOT,VSYS,IVSYS,WORK) @@ -1248,7 +1314,7 @@ SUBROUTINE EXEC(NITER,INFO,IR) & DWX*BREF/2.0, DWY*CREF/2.0, DWZ*BREF/2.0, & (DDC(K), K=1, NCONTROL) end if - 1905 FORMAT(1X,I3,20E11.3) + 1905 FORMAT(1X,I3,40E11.3) ENDIF C C------ limits on angles and rates @@ -1297,7 +1363,6 @@ SUBROUTINE EXEC(NITER,INFO,IR) C C------ set VINF() vector from new ALFA,BETA CALL VINFAB - C IF(NCONTROL.GT.0) THEN C------- set new GAM_D @@ -1343,7 +1408,7 @@ SUBROUTINE EXEC(NITER,INFO,IR) DELMAX = MAX( DELMAX , ABS(DDC(K)) ) ENDDO C - IF(DELMAX.LT.EPS) THEN + IF(DELMAX.LT.EXEC_TOL) THEN LSOL = .TRUE. C------- mark trim case as being converged ITRIM(IR) = IABS(ITRIM(IR)) @@ -1365,10 +1430,6 @@ SUBROUTINE EXEC(NITER,INFO,IR) PARVAL(IPROTZ,IR) = WROT(3)*0.5*BREF PARVAL(IPCL ,IR) = CLTOT C - -C WRITE(*,*) 'ALFA: ', ALFA -C WRITE(*,*) 'PARVAL(IPALFA,IR): ', PARVAL(IPALFA,IR) - LSEN = .TRUE. if (ltiming) then call cpu_time(t19) @@ -1380,6 +1441,7 @@ SUBROUTINE EXEC(NITER,INFO,IR) END ! EXEC + SUBROUTINE OPTGET C------------------------------------------------- C Allows toggling and setting of various @@ -1573,6 +1635,7 @@ SUBROUTINE OPTGET LAIC = .FALSE. LSRD = .FALSE. LSOL = .FALSE. + LVEL = .FALSE. ! bug 22 Mar 22 LSEN = .FALSE. C C--------------------------------- @@ -1651,8 +1714,6 @@ SUBROUTINE GETFILE(LU,FNAME) CHARACTER(*) FNAME C CHARACTER*1 ANS, DUMMY - - WRITE(*,*) LU C 1000 FORMAT(A) C @@ -1670,16 +1731,21 @@ SUBROUTINE GETFILE(LU,FNAME) WRITE(*,*) WRITE(*,*) 'File exists. Append/Overwrite/Cancel (A/O/C)? C' READ(*,1000) ANS - IF (INDEX('Aa',ANS).NE.0) THEN - 40 CONTINUE - READ(LU,1000,END=42) DUMMY - GO TO 40 - 42 CONTINUE + IF (INDEX('Aa',ANS).NE.0) THEN +C---- reopen file with append status HHY 4/17/18 + CLOSE(LU) + OPEN(LU,FILE=FNAME,STATUS='OLD',FORM='FORMATTED', + & ACCESS='APPEND',ERR=44) +C---- old append code (re-reads to EOF) +cc 40 CONTINUE +cc READ(LU,1000,END=42) DUMMY +cc GO TO 40 +cc 42 CONTINUE ELSEIF(INDEX('Oo',ANS).NE.0) THEN - REWIND(LU) + REWIND(LU) ELSE - CLOSE(LU) - LU = 0 + CLOSE(LU) + LU = 0 ENDIF RETURN C diff --git a/src/aoutmrf.f b/src/aoutmrf.f new file mode 100644 index 0000000..0709a8d --- /dev/null +++ b/src/aoutmrf.f @@ -0,0 +1,587 @@ +C*********************************************************************** +C Module: aoutmrf.f +C +C Copyright (C) 2002 Mark Drela, Harold Youngren +C +C This program is free software; you can redistribute it and/or modify +C it under the terms of the GNU General Public License as published by +C the Free Software Foundation; either version 2 of the License, or +C (at your option) any later version. +C +C This program is distributed in the hope that it will be useful, +C but WITHOUT ANY WARRANTY; without even the implied warranty of +C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +C GNU General Public License for more details. +C +C You should have received a copy of the GNU General Public License +C along with this program; if not, write to the Free Software +C Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. +C*********************************************************************** + + SUBROUTINE MRFTOT(LUN, FILEID) +C +C...PURPOSE To print out results of the vortex lattice calculation +C for the input configuration. +C +C...INPUT Configuration data for case in labeled commons +C +C...OUTPUT Machine-readable format output on logical unit LUN +C + INCLUDE 'AVL.INC' + CHARACTER*(*) FILEID + CHARACTER*50 SATYPE + + write(*,'(a,a)') 'mrftot: fileid = ', fileid +C + IF (LUN.EQ.0) RETURN +C + CALL GETSA(LNASA_SA,SATYPE,DIR) +C + CA = COS(ALFA) + SA = SIN(ALFA) +C +C---- set normalized rates in Stability or Body axes + RX_S = (WROT(1)*CA + WROT(3)*SA) * BREF/2.0 + RY_S = WROT(2) * CREF/2.0 + RZ_S = (WROT(3)*CA - WROT(1)*SA) * BREF/2.0 + RX_B = WROT(1) * BREF/2.0 + RY_B = WROT(2) * CREF/2.0 + RZ_B = WROT(3) * BREF/2.0 +C +C---- set body forces in Geometric axes + CXTOT = CDTOT*CA - CLTOT*SA + CZTOT = CDTOT*SA + CLTOT*CA +C +C---- set moments in stability axes + CRSAX = CRTOT*CA + CNTOT*SA + CMSAX = CMTOT + CNSAX = CNTOT*CA - CRTOT*SA +CCC CNSAX = CNTOT*CA - CRTOT*CA !!! Bug MD 02 Apr 04 +C +C---- dump it + WRITE(LUN,'(A)') FILEID + WRITE(LUN,'(A)') 'VERSION 1.0' +C + WRITE(LUN,'(A)') 'Vortex Lattice Output -- Total Forces' + WRITE(LUN,'(A)') TITLE(1:60) + WRITE(LUN,'(I6,6X,A)') NSURF, '| # surfaces' + WRITE(LUN,'(I6,6X,A)') NSTRIP, '| # strips' + WRITE(LUN,'(I6,6X,A)') NVOR, '| # vortices' + IF(IYSYM.GT.0) THEN + WRITE(LUN,'(ES23.15,6X,A)') YSYM, '| Y Symmetry: Wall plane' + ENDIF + IF(IYSYM.LT.0) THEN + WRITE(LUN,'(ES23.15,6X,A)') YSYM, '| Y Symmetry: Free surface' + ENDIF + IF(IZSYM.GT.0) THEN + WRITE(LUN,'(ES23.15,6X,A)') ZSYM, '| Z Symmetry: Ground plane' + ENDIF + IF(IZSYM.LT.0) THEN + WRITE(LUN,'(ES23.15,6X,A)') ZSYM, '| Z Symmetry: Free surface' + ENDIF +C + WRITE(LUN,'(3(ES23.15),6X,A)') SREF, CREF, BREF, + & '| Sref, Cref, Bref' + WRITE(LUN,'(3(ES23.15),6X,A)') XYZREF(1), XYZREF(2), XYZREF(3), + & '| Xref, Yref, Zref' +C + WRITE(LUN,'(A)') SATYPE +C + WRITE(LUN,'(A)') RTITLE(IRUN) + WRITE(LUN,'(3(ES23.15),6X,A)') ALFA/DTR, DIR*RX_B, DIR*RX_S, + & '| Alpha, pb/2V, p''b/2V' + WRITE(LUN,'(2(ES23.15),6X,A)') BETA/DTR, RY_B, + & '| Beta, qc/2V' + WRITE(LUN,'(3(ES23.15),6X,A)') AMACH , DIR*RZ_B, DIR*RZ_S, + & '| Mach, rb/2V, r''b/2V' +C + CDITOT = CDTOT - CDVTOT + WRITE(LUN,'(3(ES23.15),6X,A)') DIR*CXTOT, DIR*CRTOT, DIR*CRSAX, + & '| CXtot, Cltot, Cl''tot' + WRITE(LUN,'(2(ES23.15),6X,A)') CYTOT, CMTOT, + & '| CYtot, Cmtot' + WRITE(LUN,'(3(ES23.15),6X,A)') DIR*CZTOT, DIR*CNTOT, DIR*CNSAX, + & '| CZtot, Cntot, Cn''tot' + WRITE(LUN,'(1(ES23.15),6X,A)') CLTOT , + & '| CLtot' + WRITE(LUN,'(1(ES23.15),6X,A)') CDTOT , + & '| CDtot' + WRITE(LUN,'(2(ES23.15),6X,A)') CDVTOT, CDITOT, + & '| CDvis, CDind' + WRITE(LUN,'(4(ES23.15),6X,A)') CLFF, CDFF, CYFF, SPANEF, + & '| Trefftz Plane: CLff, CDff, CYff, e' +C + WRITE(LUN,'(A)') 'CONTROL' + WRITE(LUN,'(I6)') NCONTROL + DO K = 1, NCONTROL + WRITE(LUN,'(ES23.15,2X,A)') DELCON(K), DNAME(K) + ENDDO + + WRITE(LUN,'(A)') 'DESIGN' + WRITE(LUN,'(I6)') NDESIGN + DO K = 1, NDESIGN + WRITE(LUN,'(ES23.15,2X,A)') DELDES(K), GNAME(K) + ENDDO +C + RETURN + END ! MRFTOT + + + SUBROUTINE MRFSURF(LUN) +C +C...PURPOSE To print out surface forces from the vortex lattice calculation +C +C...INPUT Configuration data for case in labeled commons +C +C...OUTPUT Machine-readable format output for each surface on logical unit LUN +C + INCLUDE 'AVL.INC' + CHARACTER*50 SATYPE +C + IF (LUN.EQ.0) RETURN +C + CALL GETSA(LNASA_SA,SATYPE,DIR) +C +C...Print out the results +C + WRITE(LUN,'(A)') 'SURF' + WRITE(LUN,'(A)') 'VERSION 1.0' +C + WRITE(LUN,'(A)') SATYPE + WRITE(LUN,'(3(ES23.15),6X,A)') SREF, CREF, BREF, + & '| Sref, Cref, Bref' + WRITE(LUN,'(3(ES23.15),6X,A)') XYZREF(1), XYZREF(2), XYZREF(3), + & '| Xref, Yref, Zref' +C + WRITE(LUN,'(I6,6X,A)') NSURF, '| # surfaces' + DO N = 1, NSURF + CALL STRIP(STITLE(N),NT) + + WRITE(LUN,'(A)') 'SURFACE' + WRITE(LUN,'(A)') STITLE(N)(1:NT) +C +C...Force components from each surface + WRITE(LUN,'(I3,1X,9(ES23.15),1X,A,A,A)') + & N,SSURF(N),CLSURF(N),CDSURF(N),CMSURF(N), + & CYSURF(N),DIR*CNSURF(N),DIR*CRSURF(N), + & CDSURF(N)-CDVSURF(N),CDVSURF(N), + & '| Surface Forces (referred to Sref,Cref,Bref ', + & 'about Xref,Yref,Zref) : ', + & 'n Area CL CD Cm CY Cn Cl CDi CDv' +C +C--- Surface forces normalized by local reference quantities + WRITE(LUN,'(I3,1X,5(ES23.15),1X,A,A,A)') + & N,SSURF(N),CAVESURF(N), + & CL_SRF(N),CD_SRF(N), + & CDVSURF(N)*SREF/SSURF(N), + & '| Surface Forces (referred to Ssurf, Cave ', + & 'about root LE on hinge axis) : ', + & 'n Ssurf Cave cl cd cdv' + END DO + +C + RETURN + END ! MRFSURF + + + SUBROUTINE MRFBODY(LUN) +C +C...PURPOSE To print out body forces from source/doublet calculation +C +C...INPUT Configuration data for case in labeled commons +C +C...OUTPUT Machine-readable format output for each body on logical unit LUN +C + INCLUDE 'AVL.INC' + CHARACTER*50 SATYPE +C + 1000 FORMAT (A) +C + IF (LUN.EQ.0) RETURN +C + CALL GETSA(LNASA_SA,SATYPE,DIR) +C +C...Print out the results +C + WRITE(LUN,'(A)') 'BODY' + WRITE(LUN,'(A)') 'VERSION 1.0' +C + WRITE(LUN,'(A)') SATYPE + WRITE(LUN,'(3(ES23.15),6X,A)') SREF, CREF, BREF, + & '| Sref, Cref, Bref' + WRITE(LUN,'(3(ES23.15),6X,A)') XYZREF(1), XYZREF(2), XYZREF(3), + & '| Xref, Yref, Zref' +C + WRITE(LUN,'(I4,1X,A)') NBODY, '| # bodies' + DO IB = 1, NBODY + CALL STRIP(BTITLE(N),NT) + + WRITE(LUN,'(A)') 'BODY' + WRITE(LUN,'(A)') BTITLE(N)(1:NT) +C + ELBD = ELBDY(IB) + SBDY = SRFBDY(IB) + VBDY = VOLBDY(IB) + WRITE(LUN,'(I4,1X,9(ES23.15),1X,A,A,A)') + & IB,ELBD,SBDY,VBDY, + & CLBDY(IB),CDBDY(IB),CMBDY(IB), + & CYBDY(IB),DIR*CNBDY(IB),DIR*CRBDY(IB), + & '| Body Forces (referred to Sref,Cref,Bref ', + & 'about Xref,Yref,Zref) : ', + & 'Ibdy Length Asurf Vol CL CD Cm CY Cn Cl' + ENDDO +C + RETURN + END ! MRFBODY + + + SUBROUTINE MRFSTRP(LUN) +C +C...PURPOSE To print out results of the vortex lattice calculation +C for the input configuration strip and surface forces. +C +C...INPUT Configuration data for case in labeled commons +C +C...OUTPUT Machine-readable format output on logical unit LUN +C + INCLUDE 'AVL.INC' + CHARACTER*50 SATYPE +C + IF (LUN.EQ.0) RETURN +C + CALL GETSA(LNASA_SA,SATYPE,DIR) +C +C...Print out the results -> Forces by surface and strip +C + WRITE(LUN,'(A)') 'STRP' + WRITE(LUN,'(A)') 'VERSION 1.0' +C + WRITE(LUN,'(A)') SATYPE + WRITE(LUN,'(3(ES23.15),6X,A)') SREF, CREF, BREF, + & '| Sref, Cref, Bref' + WRITE(LUN,'(3(ES23.15),6X,A)') XYZREF(1), XYZREF(2), XYZREF(3), + & '| Xref, Yref, Zref' +C + WRITE(LUN,'(A,A)') 'Surface and Strip Forces by surface', + & ' (referred to Sref,Cref,Bref about Xref,Yref,Zref)' +C + WRITE(LUN,'(I6,6X,A)') NSURF, '| surfaces' + DO N = 1, NSURF + NS = NJ(N) + NV = NK(N) + J1 = JFRST(N) +C + WRITE(LUN,'(A)') 'SURFACE' + WRITE(LUN,'(A)') STITLE(N) + WRITE(LUN,'(4(I4,1X),3X,A)') + & N,NV,NS,J1, + & '| Surface #, # Chordwise, # Spanwise, First strip' + WRITE(LUN,'(2(ES23.15),3X,A)') + & SSURF(N),CAVESURF(N), + & '| Surface area Ssurf, Ave. chord Cave' +C + CDISURF = CDSURF(N)-CDVSURF(N) + WRITE(LUN,'(8(ES23.15),3X,A,A,A)') + & CLSURF(N),DIR*CRSURF(N), + & CYSURF(N), CMSURF(N), + & CDSURF(N),DIR*CNSURF(N), + & CDISURF,CDVSURF(N), + & '| CLsurf, Clsurf, CYsurf, Cmsurf, ', + & 'CDsurf, Cnsurf, CDisurf, CDvsurf', + & '; Forces referred to Sref, Cref, Bref about Xref, Yref, Zref' +C + WRITE(LUN,'(2(ES23.15),3X,A,A)') + & CL_SRF(N),CD_SRF(N), + & '| CL_srf CD_srf', + & '; Forces referred to Ssurf, Cave' +C + WRITE(LUN,'(A)') 'Strip Forces referred to Strip Area, Chord' + WRITE(LUN,'(A,A)') 'j, Xle, Yle, Zle, Chord, Area, c_cl, ai, ', + & 'cl_norm, cl, cd, cdv, cm_c/4, cm_LE, C.P.x/c' + DO JJ = 1, NS + J = J1 + JJ-1 + ASTRP = WSTRIP(J)*CHORD(J) + XCP = 999. + IF(CL_LSTRP(J).NE.0.) XCP = 0.25 - CMC4(J)/CL_LSTRP(J) + WRITE(LUN,'(I4,14(ES23.15))') + & J,RLE(1,J),RLE(2,J),RLE(3,J), + & CHORD(J),ASTRP,CNC(J),DWWAKE(J), + & CLTSTRP(J),CL_LSTRP(J),CD_LSTRP(J),CDV_LSTRP(J), + & CMC4(J),CMLE(J),XCP + END DO + END DO +C + RETURN + END ! MRFSTRP + + + SUBROUTINE MRFELE(LUN) + INCLUDE 'AVL.INC' + CHARACTER*50 SATYPE +C + IF (LUN.EQ.0) RETURN +C + CALL GETSA(LNASA_SA,SATYPE,DIR) +C + WRITE(LUN,'(A)') 'ELE' + WRITE(LUN,'(A)') 'VERSION 1.0' +C + WRITE(LUN,'(A)') SATYPE + WRITE(LUN,'(3(ES23.15),6X,A)') SREF, CREF, BREF, + & '| Sref, Cref, Bref' + WRITE(LUN,'(3(ES23.15),6X,A)') XYZREF(1), XYZREF(2), XYZREF(3), + & '| Xref, Yref, Zref' +C +C...Forces on each strip and element (long output, and slow to printout) + WRITE(LUN,'(A)') 'Vortex Strengths (by surface, by strip)' +C + WRITE(LUN,'(I6,6X,A)') NSURF, '| # surfaces' + DO N = 1, NSURF + NS = NJ(N) + NV = NK(N) + J1 = JFRST(N) +C + WRITE(LUN,'(A)') 'SURFACE' + WRITE(LUN,'(A)') STITLE(N) + WRITE(LUN,'(4(I4,1X),3X,A)') + & N,NV,NS,J1, + & '| Surface #, # Chordwise, # Spanwise, First strip' + WRITE(LUN,'(2(ES23.15),3X,A)') + & SSURF(N),CAVESURF(N), + & '| Surface area, Ave. chord' +C + CDISURF = CDSURF(N)-CDVSURF(N) + WRITE(LUN,'(8(ES23.15),3X,A,A)') + & CLSURF(N),DIR*CRSURF(N), + & CYSURF(N), CMSURF(N), + & CDSURF(N),DIR*CNSURF(N), + & CDISURF,CDVSURF(N), + & '| CLsurf, Clsurf, CYsurf, Cmsurf, ', + & 'CDsurf, Cnsurf, CDisurf, CDvsurf' + WRITE(LUN,'(2(ES23.15),3X,A,A)') + & CL_SRF(N),CD_SRF(N), + & '| CL_srf CD_srf', + & '; Forces referred to Ssurf, Cave about hinge axis thru LE' +C + DO JJ = 1, NS + J = J1 + JJ-1 + I1 = IJFRST(J) + ASTRP = WSTRIP(J)*CHORD(J) + DIHED = -ATAN2(ENSY(J),ENSZ(J))/DTR +C + WRITE(LUN,'(A)') 'STRIP' + WRITE(LUN,'(3(I4),2X,A)') J,NV,I1, + & '| Strip #, # Chordwise, First Vortex' + WRITE(LUN,'(8(ES23.15),2X,A,A)') + & RLE(1,J),CHORD(J),AINC(J)/DTR, + & RLE(2,J),WSTRIP(J),ASTRP, + & RLE(3,J),DIHED, + & '| Xle, Ave. Chord, Incidence (deg), Yle, ', + & 'Strip Width, Strip Area, Zle, Strip Dihed (deg)' + WRITE(LUN,'(9(ES23.15),2X,A)') + & CL_LSTRP(J), CD_LSTRP(J), CDV_LSTRP(J), + & CNRMSTRP(J), CAXLSTRP(J), + & CNC(J), DWWAKE(J), + & CMLE(J), CMC4(J), + & '| cl, cd, cdv, cn, ca, cnc, wake dnwsh, cmLE, cm c/4' +C + !!WRITE(LUN,'(I4,2X,A)') NV, '| # vortices' + DO II = 1, NV + I = I1 + (II-1) + XM = 0.5*(RV1(1,I)+RV2(1,I)) + YM = 0.5*(RV1(2,I)+RV2(2,I)) + ZM = 0.5*(RV1(3,I)+RV2(3,I)) + WRITE(LUN,'(I4,6(ES23.15),2X,A)') + & I,XM,YM,ZM,DXV(I),SLOPEC(I),DCP(I), + & '| I, X, Y, Z, DX, Slope, dCp' + END DO + END DO + END DO +C + END ! MRFELE + + + SUBROUTINE MRFHINGE(LUN) +C +C...PURPOSE To print out results of the vortex lattice calculation +C for the input configuration. +C +C...INPUT Configuration data for case in labeled commons +C +C...OUTPUT Machine-readable format output on logical unit LUN +C + INCLUDE 'AVL.INC' + CHARACTER*50 SATYPE +C + IF (LUN.EQ.0) RETURN +C + CALL GETSA(LNASA_SA,SATYPE,DIR) +C +C...Print out the results +C + WRITE(LUN,'(A)') 'HINGE' + WRITE(LUN,'(A)') 'VERSION 1.0' +C + WRITE(LUN,'(A)') SATYPE + WRITE(LUN,'(2(ES23.15),6X,A)') SREF, CREF, + & '| Sref, Cref' +C +C...Hinge moments for each CONTROL + !!WRITE(LUN,'(A)') 'Control Hinge Moments (referred to Sref, Cref)' + WRITE(LUN,'(I4,2X,A)') NCONTROL, '| # controls' + DO N = 1, NCONTROL + WRITE(LUN,'(ES23.15,2X,A,2X,A,A)') CHINGE(N), DNAME(N), + & '| Control Hinge Moments (referred to Sref, Cref) : ', + & 'Chinge, Control' + END DO +C + RETURN + END ! MRFHINGE + + + SUBROUTINE MRFCNC(LUN) +C +C...PURPOSE To write out a CNC loading file +C for the input configuration strips +C +C...INPUT Configuration data for case in labeled commons +C +C...OUTPUT Machine-readable format output on logical unit LUN +C + INCLUDE 'AVL.INC' + CHARACTER*1 ANS + CHARACTER*256 FNAM + SAVE FNAM + DATA FNAM /' '/ +C + IF (LUN.EQ.0) RETURN +C +C...Print out the results -> strip loadings +C + WRITE(LUN,'(A)') 'CNC' + WRITE(LUN,'(A)') 'VERSION 1.0' +C + WRITE(LUN,'(A,A)') 'Strip Loadings: ', + & ' XM, YM, ZM, CNCM, CLM, CHM, DYM, ASM' + WRITE(LUN,'(I4,3X,A)') NSTRIP, '| # strips' +C + DO J=1, NSTRIP + I = IJFRST(J) + XM = 0.5*(RV1(1,I)+RV2(1,I)) + YM = 0.5*(RV1(2,I)+RV2(2,I)) + ZM = 0.5*(RV1(3,I)+RV2(3,I)) + CNCM = CNC(J) + CLM = CL_LSTRP(J) + CHM = CHORD(J) + DYM = WSTRIP(J) + ASM = DYM*CHM + WRITE(LUN,'(8(ES23.15))') XM,YM,ZM,CNCM,CLM,CHM,DYM,ASM + END DO +C + RETURN + END ! MRFCNC + + + SUBROUTINE MRFVM(LU) +C...PURPOSE Print STRIP SHEAR and BENDING QUANTITIES, ie. V, BM +C Integrates spanload to get shear and bending moment +C NOTE: only works for single surface at at time (ie. V,BM on each panel) +C + INCLUDE 'AVL.INC' + REAL V(NSMAX), BM(NSMAX), YSTRP(NSMAX) + CHARACTER*256 FNAMVB +C + IF (LU.EQ.0) RETURN +C + WRITE(LU,'(A)') 'VM' + WRITE(LU,'(A)') 'VERSION 1.0' +C + WRITE(LU,'(A)') 'Shear/q and Bending Moment/q vs Y' + WRITE(LU,'(A)') TITLE(1:60) + WRITE(LU,'(6(ES23.15),2X,A)') + & AMACH,ALFA/DTR,CLTOT,BETA/DTR,SREF,BREF, + & '| Mach, alpha, CLtot, beta, Sref, Bref' +C +C---- Process the surfaces one by one, calculating shear and bending on each, +C with moments refered to Y=0 (centerline) +C + WRITE(LU,'(I6,6X,A)') NSURF, '| # surfaces' +C + DO N = 1, NSURF + J1 = JFRST(N) + JN = J1 + NJ(N) - 1 +C + WRITE(LU,'(A)') 'SURFACE' + WRITE(LU,'(A)') STITLE(N) + WRITE(LU,'(2(I4,1X),3X,A)') + & N, NJ(N), + & '| Surface #, # strips' +C + YMIN = 1.0E10 + YMAX = -1.0E10 + DO J = J1, JN + IV = IJFRST(J) + YMIN = MIN(YMIN,RV1(2,IV),RV2(2,IV)) + YMAX = MAX(YMAX,RV1(2,IV),RV2(2,IV)) + ENDDO +C +C------ Integrate spanload from first strip to last strip defined for +C this surface to get shear and bending moment + CNCLST = 0.0 + BMLST = 0.0 + WLST = 0.0 + VLST = 0.0 +C + JF = J1 + JL = JN + JINC = 1 +C +C------ Integrate from first to last strip in surface + DO J = JL, JF, -JINC + JJ = JINC*(J - JF + JINC) +C + DY = 0.5*(WSTRIP(J)+WLST) + YSTRP(JJ) = RLE(2,J) + V(JJ) = VLST + 0.5*(CNC(J)+CNCLST) * DY + BM(JJ) = BMLST + 0.5*(V(JJ)+VLST) * DY +C + VLST = V(JJ) + BMLST = BM(JJ) + CNCLST = CNC(J) + WLST = WSTRIP(J) + ENDDO +C +C------ Inboard edge Y,Vz,Mx + VROOT = VLST + CNCLST * 0.5*DY + BMROOT = BMLST + 0.5*(VROOT+VLST) * 0.5*DY + VTIP = 0.0 + BMTIP = 0.0 + IF(IMAGS(N).GE.0) THEN + YROOT = RLE1(2,J1) + YTIP = RLE2(2,JN) + ELSE + YROOT = RLE2(2,J1) + YTIP = RLE1(2,JN) + ENDIF +C + DIR = 1.0 + IF(YMIN+YMAX.LT.0.0) DIR = -1.0 +C + WRITE(LU,'(2(ES23.15),2X,A)') 2.0*YMIN/BREF, 2.0*YMAX/BREF, + & '| 2Ymin/Bref, 2Ymax/Bref' +C + WRITE(LU,'(3(ES23.15),2X,A)') + & 2.*YROOT/BREF,VROOT/SREF,DIR*BMROOT/SREF/BREF, + & '| 2Y/Bref, Vz/(q*Sref), Mx/(q*Bref*Sref) : root' + DO J = 1, NJ(N) + WRITE(LU,'(3(ES23.15),2X,A)') + & 2.*YSTRP(J)/BREF,V(J)/SREF,DIR*BM(J)/SREF/BREF, + & '| 2Y/Bref, Vz/(q*Sref), Mx/(q*Bref*Sref)' + ENDDO + WRITE(LU,'(3(ES23.15),2X,A)') + & 2.*YTIP/BREF,VTIP/SREF,DIR*BMTIP/SREF/BREF, + & '| 2Y/Bref, Vz/(q*Sref), Mx/(q*Bref*Sref) : tip' + ENDDO +C + RETURN + END ! MRFVM diff --git a/src/aoutput.f b/src/aoutput.f index 0b41915..7a3878f 100644 --- a/src/aoutput.f +++ b/src/aoutput.f @@ -5,7 +5,7 @@ C C This program is free software; you can redistribute it and/or modify C it under the terms of the GNU General Public License as published by -C the Free Software Foundation; either version 2 of the License, or +C the Free Software Foundation; either version 2 of the License, or C (at your option) any later version. C C This program is distributed in the hope that it will be useful, @@ -127,8 +127,8 @@ SUBROUTINE OUTTOT(LUN) & /2X,'CZtot =',F10.5,5X,'Cntot =',F10.5,5X,'Cn''tot =',F10.5 &//2X,'CLtot =',F10.5 & /2X,'CDtot =',F10.5 - & /2X,'CDvis =',F10.5,5X,'CDind =',F10.5 - & /2X,'CLff =',F10.5,5X,'CDff =',F10.5,4X,'| Trefftz' + & /2X,'CDvis =',F10.5,5X,'CDind =',F10.7 + & /2X,'CLff =',F10.5,5X,'CDff =',F10.7,4X,'| Trefftz' & /2X,'CYff =',F10.5,5X,' e =',F10.4,4X,'| Plane ' ) C 231 FORMAT(3X,A,'=',F10.5) @@ -180,7 +180,7 @@ SUBROUTINE OUTSURF(LUN) CALL STRIP(STITLE(N),NT) WRITE (LUN,221) N,SSURF(N),CAVESURF(N), & CL_SRF(N),CD_SRF(N), - & CDVSURF(N)*SREF/SSURF(N),CMLE_SRF(N), + & CDVSURF(N)*SREF/SSURF(N), & STITLE(N)(1:NT) END DO WRITE(LUN,200) @@ -198,8 +198,8 @@ SUBROUTINE OUTSURF(LUN) 220 FORMAT (/' Surface Forces (referred to Ssurf, Cave ', & 'about root LE on hinge axis)'// & 2X,' n',5X,'Ssurf',6X,'Cave', - & 7X,'cl',7X,'cd',6X,'cdv',4x,'cm_LE') - 221 FORMAT (2X,I2,F10.3,F10.3,4(1X,F8.4),2X,A) + & 7X,'cl',7X,'cd',6X,'cdv') + 221 FORMAT (2X,I2,F10.3,F10.3,3(1X,F8.4),2X,A) C RETURN END ! OUTSURF @@ -249,10 +249,14 @@ SUBROUTINE OUTBODY(LUN) & /' ',A // & 5X,'Sref =',G12.4, 3X,'Cref =',F10.4,3X,'Bref =',F10.4/ & 5X,'Xref =',2X,F10.4,3X,'Yref =',F10.4,3X,'Zref =',F10.4// - & 'Ibdy',4X,'Length',5X,'Asurf',7X,'Vol', - & 6X,'CL',6X,'CD',6X,'Cm', - & 6X,'CY',6X,'Cn',6X,'Cl') - 211 FORMAT (I4,3(1X,F9.3),6F8.4,3X,A) +c & 'Ibdy',4X,'Length',5X,'Asurf',7X,'Vol', +c & 6X,'CL',6X,'CD',6X,'Cm', +c & 6X,'CY',6X,'Cn',6X,'Cl') +c 211 FORMAT (I4,3(1X,F9.3),6F8.4,3X,A) + & 'Ibdy',7X,'Length',8X,'Asurf',10X,'Vol', + & 10X,'CL',10X,'CD',10X,'Cm', + & 10X,'CY',10X,'Cn',10X,'Cl') + 211 FORMAT (I4,3(1X,F12.6),6F12.6,3X,A) 212 FORMAT (/) C RETURN @@ -280,8 +284,9 @@ SUBROUTINE OUTSTRP(LUN) C C...Print out the results -> Forces by surface and strip WRITE(LUN,200) - WRITE(LUN,210) - WRITE(LUN,212) SATYPE + WRITE(LUN,205) + WRITE (LUN,206) SREF,CREF,BREF, + & XYZREF(1), XYZREF(2), XYZREF(3) DO N = 1, NSURF NS = NJ(N) NV = NK(N) @@ -289,6 +294,7 @@ SUBROUTINE OUTSTRP(LUN) C WRITE (LUN,211) N,STITLE(N),NV,NS,J1,SSURF(N),CAVESURF(N) CDISURF = CDSURF(N)-CDVSURF(N) + WRITE (LUN,212) SATYPE WRITE (LUN,213) CLSURF(N),DIR*CRSURF(N), & CYSURF(N), CMSURF(N), & CDSURF(N),DIR*CNSURF(N), @@ -299,46 +305,48 @@ SUBROUTINE OUTSTRP(LUN) J = J1 + JJ-1 ASTRP = WSTRIP(J)*CHORD(J) XCP = 999. -cc IF(CLSTRP(J).NE.0.) XCP = 0.25 - CMC4(J)/CLSTRP(J) !!! BUG 13 Jan 12 MD IF(CL_LSTRP(J).NE.0.) XCP = 0.25 - CMC4(J)/CL_LSTRP(J) IF(XCP.LT.2.0 .AND. XCP.GT.-1.0) THEN WRITE (LUN,217) - & J,RLE(2,J),CHORD(J),ASTRP,CNC(J),DWWAKE(J), + & J,RLE(1,J),RLE(2,J),RLE(3,J), + & CHORD(J),ASTRP,CNC(J),DWWAKE(J), & CLTSTRP(J), CL_LSTRP(J),CD_LSTRP(J),CDV_LSTRP(J), & CMC4(J),CMLE(J),XCP ELSE WRITE (LUN,217) - & J,RLE(2,J),CHORD(J),ASTRP,CNC(J),DWWAKE(J), + & J,RLE(1,J),RLE(2,J),RLE(3,J), + & CHORD(J),ASTRP,CNC(J),DWWAKE(J), & CLTSTRP(J),CL_LSTRP(J),CD_LSTRP(J),CDV_LSTRP(J), - & CMC4(J),CMLE(J) + & CMC4(J),CMLE(J),999. ENDIF END DO END DO - WRITE(LUN,200) + WRITE (LUN,200) C 200 FORMAT(1X, &'---------------------------------------------------------------') - 210 FORMAT (' Surface and Strip Forces by surface') + 205 FORMAT (' Surface and Strip Forces by surface') + 206 FORMAT(/2X, 'Sref =',G12.5,3X,'Cref =',G12.5,3X,'Bref =',G12.5 + & /2X, 'Xref =',G12.5,3X,'Yref =',G12.5,3X,'Zref =',G12.5 ) +C 211 FORMAT (/2X,'Surface #',I2,5X,A/ & 5X,'# Chordwise =',I3,3X,'# Spanwise =',I3, & 5X,'First strip =',I3/ - & 5X,'Surface area =',F12.6,5X,' Ave. chord =',F12.6) - 212 FORMAT (/' Forces referred to Sref, Cref, Bref ', - & 'about Xref, Yref, Zref'/ - & ' ',A) + & 5X,'Surface area Ssurf =',F12.6, + & 5X,'Ave. chord Cave =',F12.6) + 212 FORMAT (/' Forces referred to Sref, Cref, Bref ', + & 'about Xref, Yref, Zref'/' ',A) 213 FORMAT ( 5X,'CLsurf =',F10.5,5X,'Clsurf =',F10.5, & /5X,'CYsurf =',F10.5,5X,'Cmsurf =',F10.5, & /5X,'CDsurf =',F10.5,5X,'Cnsurf =',F10.5, & /5X,'CDisurf =',F10.5,5x,'CDvsurf =',F10.5) - 214 FORMAT (/' Forces referred to Ssurf, Cave ', - & 'about hinge axis thru LE'/ - & 5X,'CLsurf =',F10.5,5X,'CDsurf =',F10.5/ - & 5X,'Deflect =',F10.5,5X,'CmLEsurf=',F10.5) + 214 FORMAT (/' Forces referred to Ssurf, Cave '/ + & 5X,'CLsurf =',F10.5,5X,'CDsurf =',F10.5) 216 FORMAT (/' Strip Forces referred to Strip Area, Chord'/ - & 2X,' j ',5X,'Yle',4X,'Chord',5X,'Area', - & 5X,'c cl',6X,'ai',6X,'cl_norm',2X,'cl',7X,'cd',7X, - & 'cdv',4x,'cm_c/4',4x,'cm_LE',2x,'C.P.x/c') - 217 FORMAT (2X,I4,11(1X,F8.4),1X,F8.3) + & 2X,' j ',4X,'Xle',6X,'Yle',6X,'Zle',6X,'Chord',4X,'Area', + & 5X,'c_cl',5X,'ai',5X,'cl_norm',4X,'cl',7X,'cd',7X, + & 'cdv',4x,'cm_c/4',5x,'cm_LE',3x,'C.P.x/c') + 217 FORMAT (2X,I4,13(1X,F8.4),1X,F8.3) C RETURN END ! OUTSTRP @@ -404,8 +412,8 @@ SUBROUTINE OUTELE(LUN) 230 FORMAT (' Vortex Strengths (by surface, by strip)') 231 FORMAT (/1X,78('*')/2X,'Surface #',I2,5X,A/ & 5X,'# Chordwise =',I3,3X,'# Spanwise =',I3, - & 3X,'First strip =',I3/ - & 5X,'Surface area =',F12.6,5X,' Ave. chord =',F12.6) + & 3X,'First strip =',I4/ + & 4X,'Surface area =',F12.6,5X,' Ave. chord =',F12.6) C 212 FORMAT (/' Forces referred to Sref, Cref, Bref ', & 'about Xref, Yref, Zref'/ @@ -432,7 +440,7 @@ SUBROUTINE OUTELE(LUN) & /4X,'cmLE=',F10.5,4X,'cm c/4 =',F10.5, & //4X,'I',8X,'X ',8X,'Y ',8X,'Z ',8X,'DX ', & 6X,'Slope',8X,'dCp') - 234 FORMAT (2X,I3,6(2X,F10.5)) + 234 FORMAT (1X,I4,6(2X,F10.5)) C END ! OUTELE @@ -492,7 +500,7 @@ SUBROUTINE OUTCNC(LUN) C INCLUDE 'AVL.INC' CHARACTER*1 ANS - CHARACTER*120 FNAM + CHARACTER*256 FNAM SAVE FNAM DATA FNAM /' '/ C @@ -526,7 +534,8 @@ SUBROUTINE OUTCNC(LUN) END DO C 210 FORMAT (//' *** Strip Loadings') - 212 FORMAT (3(F8.3,1X),2(F10.4,1X),2(F8.4,1X),F9.4) +ccc 212 FORMAT (3(F8.3,1X),2(F10.4,1X),2(F8.4,1X),F9.4) + 212 FORMAT (8(F12.6,1X)) 213 FORMAT (F8.3,1X,3(F10.4,1X)) C RETURN @@ -534,12 +543,15 @@ SUBROUTINE OUTCNC(LUN) - SUBROUTINE DERMATM(LU) + SUBROUTINE DERMATM(LU, USEMRF) C--------------------------------------------------------- C Calculates and outputs stability derivative matrix C for current ALFA, BETA. +C Forces and rotations are about stability axes, +C moments are about body axes C--------------------------------------------------------- INCLUDE 'AVL.INC' + LOGICAL USEMRF CHARACTER*50 SATYPE REAL WROT_RX(3), WROT_RZ(3), WROT_A(3) C @@ -575,7 +587,6 @@ SUBROUTINE DERMATM(LU) WROT_A(2) = 0. WROT_A(3) = -RZ*SA + RX*CA !!! = WROT(1) C -C C---- set force derivatives in stability axes CL_AL = CLTOT_U(1)*VINF_A(1) + CLTOT_U(4)*WROT_A(1) & + CLTOT_U(2)*VINF_A(2) + CLTOT_U(5)*WROT_A(2) @@ -627,7 +638,132 @@ SUBROUTINE DERMATM(LU) CN_RY = CNTOT_U(5) CN_RZ = CNTOT_U(6)*WROT_RZ(3) + CNTOT_U(4)*WROT_RZ(1) C +C... machine-readable output format + IF(USEMRF) THEN + CALL MRFTOT(LU, 'DERMATM') +C + WRITE(LU,'(/A)') 'Derivatives...' +C + WRITE(LU,'(A)') 'alpha, beta' + WRITE(LU,'(2(ES23.15),2X,A)') + & CL_AL, CL_BE, + & '| z force CL : CLa, CLb' + WRITE(LU,'(2(ES23.15),2X,A)') + & CY_AL, CY_BE, + & '| y force CY : CYa, CYb' + WRITE(LU,'(2(ES23.15),2X,A)') + & DIR*CR_AL, DIR*CR_BE, + & '| roll x mom. : Cla, Clb' + WRITE(LU,'(2(ES23.15),2X,A)') + & CM_AL, CM_BE, + & '| pitch y mom. : Cma, Cmb' + WRITE(LU,'(2(ES23.15),2X,A)') + & DIR*CN_AL, DIR*CN_BE, + & '| yaw z mom. : Cna, Cnb' +C + WRITE(LU,'(A)') 'roll rate p, pitch rate q, yaw rate r' + WRITE(LU,'(3(ES23.15),2X,A)') + & CL_RX*2.0/BREF, + & CL_RY*2.0/CREF, + & CL_RZ*2.0/BREF, + & '| z force : CLp, CLq, CLr' + WRITE(LU,'(3(ES23.15),2X,A)') + & CY_RX*2.0/BREF, + & CY_RY*2.0/CREF, + & CY_RZ*2.0/BREF, + & '| y force : CYp, CYq, CYr' + WRITE(LU,'(3(ES23.15),2X,A)') + & DIR*CR_RX*2.0/BREF, + & DIR*CR_RY*2.0/CREF, + & DIR*CR_RZ*2.0/BREF, + & '| roll x mom. : Clp, Clq, Clr' + WRITE(LU,'(3(ES23.15),2X,A)') + & CM_RX*2.0/BREF, + & CM_RY*2.0/CREF, + & CM_RZ*2.0/BREF, + & '| pitch y mom. : Cmp, Cmq, Cmr' + WRITE(LU,'(3(ES23.15),2X,A)') + & DIR*CN_RX*2.0/BREF, + & DIR*CN_RY*2.0/CREF, + & DIR*CN_RZ*2.0/BREF, + & '| yaw z mom. : Cnp, Cnq, Cnr' +C + WRITE(LU,'(I4,2X,A)') NCONTROL, '| # control vars' + IF(NCONTROL.GT.0) THEN + DO K=1, NCONTROL + WRITE(LU,'(A)') DNAME(K) + ENDDO + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & ( CLTOT_D(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| z force : CLd*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & ( CYTOT_D(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| y force : CYd*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & (DIR*CRTOT_D(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| roll x mom. : Cld*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & ( CMTOT_D(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| pitch y mom. : Cmd*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & (DIR*CNTOT_D(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| yaw z mom. : Cnd*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & ( CDFF_D(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| Trefftz drag : CDffd*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & ( SPANEF_D(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| span eff. : ed*' + ENDIF +C + WRITE(LU,'(I4,2X,A)') NDESIGN, '| # design vars' + IF(NDESIGN.GT.0) THEN + DO K=1, NDESIGN + WRITE(LU,'(A)') GNAME(K) + ENDDO + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & ( CLTOT_G(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| z force : CLg*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & ( CYTOT_G(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| y force : CYg*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & (DIR*CRTOT_G(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| roll x mom. : Clg*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & ( CMTOT_G(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| pitch y mom. : Cmg*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & (DIR*CNTOT_G(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| yaw z mom. : Cng*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & ( CDFF_G(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| Trefftz drag : CDffg*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & ( SPANEF_G(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| span eff. : eg*' + ENDIF C + IF(CL_AL .NE. 0.0) THEN + XNP = XYZREF(1) - CREF*CM_AL/CL_AL + ELSE + XNP = -1.0E30 + ENDIF + WRITE(LU,'(ES23.15,2X,A)') XNP, + & '| Neutral point Xnp' +C + IF(ABS(CR_RZ*CN_BE) .GT. 0.0001) THEN + BB = CR_BE*CN_RZ / (CR_RZ*CN_BE) + ELSE + BB = -1.0E30 + ENDIF + WRITE(LU,'(ES23.15,2X,A)') BB, + & '| Clb Cnr / Clr Cnb ( > 1 if spirally stable )' +C + RETURN + ENDIF +C +C... original output CALL OUTTOT(LU) C WRITE(LU,7004) @@ -701,30 +837,30 @@ SUBROUTINE DERMATM(LU) IF(NCONTROL.GT.0) THEN C WRITE(LU,8106) (DNAME(K), K, K=1, NCONTROL) - 8106 FORMAT(/14X,20(4X,A12, ' d',I1,' ')) + 8106 FORMAT(/14X,20(4X,A12, ' d',I2.2,' ')) WRITE(LU,8107) (' ',K=1, NCONTROL) 8107 FORMAT( 14X,20(3X,A,'----------------')) C WRITE(LU,8110) (' ',K,CLTOT_D(K), K=1, NCONTROL) - 8110 FORMAT(' z force |',20(A,' CLd',I1,' =',F11.6)) + 8110 FORMAT(' z force |',20(A,' CLd',I2.2,' =',F11.6)) C WRITE(LU,8120) (' ',K,CYTOT_D(K), K=1, NCONTROL) - 8120 FORMAT(' y force |',20(A,' CYd',I1,' =',F11.6)) + 8120 FORMAT(' y force |',20(A,' CYd',I2.2,' =',F11.6)) C - WRITE(LU,8140) (' ',K,CRTOT_D(K), K=1, NCONTROL) - 8140 FORMAT(' roll x mom.|',20(A,' Cld',I1,' =',F11.6)) + WRITE(LU,8140) (' ',K,DIR*CRTOT_D(K), K=1, NCONTROL) + 8140 FORMAT(' roll x mom.|',20(A,' Cld',I2.2,' =',F11.6)) C WRITE(LU,8150) (' ',K, CMTOT_D(K), K=1, NCONTROL) - 8150 FORMAT(' pitch y mom.|',20(A,' Cmd',I1,' =',F11.6)) + 8150 FORMAT(' pitch y mom.|',20(A,' Cmd',I2.2,' =',F11.6)) C - WRITE(LU,8160) (' ',K,CNTOT_D(K), K=1, NCONTROL) - 8160 FORMAT(' yaw z mom.|',20(A,' Cnd',I1,' =',F11.6)) + WRITE(LU,8160) (' ',K,DIR*CNTOT_D(K), K=1, NCONTROL) + 8160 FORMAT(' yaw z mom.|',20(A,' Cnd',I2.2,' =',F11.6)) C WRITE(LU,8170) (' ',K, CDFF_D(K), K=1, NCONTROL) - 8170 FORMAT(' Trefftz drag|',20(A,'CDffd',I1,' =',F11.6)) + 8170 FORMAT(' Trefftz drag|',20(A,'CDffd',I2.2,' =',F11.6)) C WRITE(LU,8180) (' ',K, SPANEF_D(K), K=1, NCONTROL) - 8180 FORMAT(' span eff. |',20(A,' ed',I1,' =',F11.6)) + 8180 FORMAT(' span eff. |',20(A,' ed',I2.2,' =',F11.6)) C WRITE(LU,*) WRITE(LU,*) @@ -734,30 +870,30 @@ SUBROUTINE DERMATM(LU) IF(NDESIGN.GT.0) THEN C WRITE(LU,8206) (GNAME(K), K, K=1, NDESIGN) - 8206 FORMAT(/14X,20(4X,A12, ' g',I1,' ')) + 8206 FORMAT(/14X,20(4X,A12, ' g',I2.2,' ')) WRITE(LU,8207) (' ',K=1, NDESIGN) 8207 FORMAT( 14X,20(3X,A,'----------------')) C WRITE(LU,8210) (' ',K,CLTOT_G(K), K=1, NDESIGN) - 8210 FORMAT(' z force |',20(A,' CLg',I1,' =',F11.6)) + 8210 FORMAT(' z force |',20(A,' CLg',I2.2,' =',F11.6)) C WRITE(LU,8220) (' ',K,CYTOT_G(K), K=1, NDESIGN) - 8220 FORMAT(' y force |',20(A,' CYg',I1,' =',F11.6)) + 8220 FORMAT(' y force |',20(A,' CYg',I2.2,' =',F11.6)) C WRITE(LU,8230) (' ',K,DIR*CRTOT_G(K), K=1, NDESIGN) - 8230 FORMAT(' roll x mom.|',20(A,' Clg',I1,' =',F11.6)) + 8230 FORMAT(' roll x mom.|',20(A,' Clg',I2.2,' =',F11.6)) C WRITE(LU,8240) (' ',K, CMTOT_G(K), K=1, NDESIGN) - 8240 FORMAT(' pitch y mom.|',20(A,' Cmg',I1,' =',F11.6)) + 8240 FORMAT(' pitch y mom.|',20(A,' Cmg',I2.2,' =',F11.6)) C WRITE(LU,8250) (' ',K,DIR*CNTOT_G(K), K=1, NDESIGN) - 8250 FORMAT(' yaw z mom.|',20(A,' Cng',I1,' =',F11.6)) + 8250 FORMAT(' yaw z mom.|',20(A,' Cng',I2.2,' =',F11.6)) C WRITE(LU,8260) (' ',K, CDFF_G(K), K=1, NDESIGN) - 8260 FORMAT(' Trefftz drag|',20(A,'CDffg',I1,' =',F11.6)) + 8260 FORMAT(' Trefftz drag|',20(A,'CDffg',I2.2,' =',F11.6)) C WRITE(LU,8270) (' ',K, SPANEF_G(K), K=1, NDESIGN) - 8270 FORMAT(' span eff. |',20(A,' eg',I1,' =',F11.6)) + 8270 FORMAT(' span eff. |',20(A,' eg',I2.2,' =',F11.6)) C WRITE(LU,*) WRITE(LU,*) @@ -776,19 +912,20 @@ SUBROUTINE DERMATM(LU) 8402 FORMAT(/' Clb Cnr / Clr Cnb =', F11.6, & ' ( > 1 if spirally stable )') ENDIF - C RETURN END ! DERMATM - - SUBROUTINE DERMATS(LU) + SUBROUTINE DERMATS(LU, USEMRF) C--------------------------------------------------------- C Calculates and outputs stability derivative matrix C for current ALFA, BETA. +C Forces and rotations are about stability axes, +C moments are about stability axes C--------------------------------------------------------- INCLUDE 'AVL.INC' + LOGICAL USEMRF CHARACTER*50 SATYPE REAL WROT_RX(3), WROT_RZ(3), WROT_A(3) REAL @@ -904,7 +1041,132 @@ SUBROUTINE DERMATS(LU) CN_RY = CNSAX_U(5) CN_RZ = CNSAX_U(6)*WROT_RZ(3) + CNSAX_U(4)*WROT_RZ(1) C +C... machine-readable output format + IF(USEMRF) THEN + CALL MRFTOT(LU, 'DERMATS') +C + WRITE(LU,'(/A)') 'Stability-axis derivatives...' +C + WRITE(LU,'(A)') 'alpha, beta' + WRITE(LU,'(2(ES23.15),2X,A)') + & CL_AL, CL_BE, + & '| z'' force CL : CLa, CLb' + WRITE(LU,'(2(ES23.15),2X,A)') + & CY_AL, CY_BE, + & '| y force CY : CYa, CYb' + WRITE(LU,'(2(ES23.15),2X,A)') + & DIR*CR_AL, DIR*CR_BE, + & '| x'' mom. Cl'' : Cla, Clb' + WRITE(LU,'(2(ES23.15),2X,A)') + & CM_AL, CM_BE, + & '| y mom. Cm : Cma, Cmb' + WRITE(LU,'(2(ES23.15),2X,A)') + & DIR*CN_AL, DIR*CN_BE, + & '| z'' mom. Cn'' : Cna, Cnb' +C + WRITE(LU,'(A)') 'roll rate p'', pitch rate q'', yaw rate r''' + WRITE(LU,'(3(ES23.15),2X,A)') + & CL_RX*2.0/BREF, + & CL_RY*2.0/CREF, + & CL_RZ*2.0/BREF, + & '| z'' force CL : CLp, CLq, CLr' + WRITE(LU,'(3(ES23.15),2X,A)') + & CY_RX*2.0/BREF, + & CY_RY*2.0/CREF, + & CY_RZ*2.0/BREF, + & '| y force CY : CYp, CYq, CYr' + WRITE(LU,'(3(ES23.15),2X,A)') + & DIR*CR_RX*2.0/BREF, + & DIR*CR_RY*2.0/CREF, + & DIR*CR_RZ*2.0/BREF, + & '| x'' mom. Cl'' : Clp, Clq, Clr' + WRITE(LU,'(3(ES23.15),2X,A)') + & CM_RX*2.0/BREF, + & CM_RY*2.0/CREF, + & CM_RZ*2.0/BREF, + & '| y mom. Cm : Cmp, Cmq, Cmr' + WRITE(LU,'(3(ES23.15),2X,A)') + & DIR*CN_RX*2.0/BREF, + & DIR*CN_RY*2.0/CREF, + & DIR*CN_RZ*2.0/BREF, + & '| z'' mom. Cn'' : Cnp, Cnq, Cnr' +C + WRITE(LU,'(I4,2X,A)') NCONTROL, '| # control vars' + IF(NCONTROL.GT.0) THEN + DO K=1, NCONTROL + WRITE(LU,'(A)') DNAME(K) + ENDDO + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & ( CLTOT_D(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| z'' force CL : CLd*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & ( CYTOT_D(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| y force CY : CYd*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & (DIR*CRTOT_D(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| x'' mom. Cl'' : Cld*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & ( CMTOT_D(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| y mom. Cm : Cmd*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & (DIR*CNTOT_D(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| z'' mom. Cn : Cnd*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & ( CDFF_D(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| Trefftz drag : CDffd*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & ( SPANEF_D(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| span eff. : ed*' + ENDIF +C + WRITE(LU,'(I4,2X,A)') NDESIGN, '| # design vars' + IF(NDESIGN.GT.0) THEN + DO K=1, NDESIGN + WRITE(LU,'(A)') GNAME(K) + ENDDO + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & ( CLTOT_G(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| z'' force CL : CLg*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & ( CYTOT_G(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| y force CY : CYg*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & (DIR*CRTOT_G(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| x'' mom. Cl'' : Clg*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & ( CMTOT_G(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| y mom. Cm : Cmg*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & (DIR*CNTOT_G(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| z'' mom. Cn : Cng*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & ( CDFF_G(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| Trefftz drag : CDffg*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & ( SPANEF_G(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| span eff. : eg*' + ENDIF +C + IF(CL_AL .NE. 0.0) THEN + XNP = XYZREF(1) - CREF*CM_AL/CL_AL + ELSE + XNP = -1.0E30 + ENDIF + WRITE(LU,'(ES23.15,2X,A)') XNP, + & '| Neutral point Xnp' +C + IF(ABS(CR_RZ*CN_BE) .GT. 0.0001) THEN + BB = CR_BE*CN_RZ / (CR_RZ*CN_BE) + ELSE + BB = -1.0E30 + ENDIF + WRITE(LU,'(ES23.15,2X,A)') BB, + & '| Clb Cnr / Clr Cnb ( > 1 if spirally stable )' +C + RETURN + ENDIF C +C... original output CALL OUTTOT(LU) C WRITE(LU,7004) @@ -978,30 +1240,30 @@ SUBROUTINE DERMATS(LU) IF(NCONTROL.GT.0) THEN C WRITE(LU,8106) (DNAME(K), K, K=1, NCONTROL) - 8106 FORMAT(/14X,20(4X,A12, ' d',I1,' ')) + 8106 FORMAT(/14X,20(4X,A12, ' d',I2.2,' ')) WRITE(LU,8107) (' ',K=1, NCONTROL) 8107 FORMAT( 14X,20(3X,A,'----------------')) C WRITE(LU,8110) (' ',K,CLTOT_D(K), K=1, NCONTROL) - 8110 FORMAT(' z'' force CL |' ,20(A,' CLd',I1,' =',F11.6)) + 8110 FORMAT(' z'' force CL |' ,20(A,' CLd',I2.2,' =',F11.6)) C WRITE(LU,8120) (' ',K,CYTOT_D(K), K=1, NCONTROL) - 8120 FORMAT(' y force CY |' ,20(A,' CYd',I1,' =',F11.6)) + 8120 FORMAT(' y force CY |' ,20(A,' CYd',I2.2,' =',F11.6)) C - WRITE(LU,8140) (' ',K,CRTOT_D(K), K=1, NCONTROL) - 8140 FORMAT(' x'' mom. Cl''|',20(A,' Cld',I1,' =',F11.6)) + WRITE(LU,8140) (' ',K,DIR*CRTOT_D(K), K=1, NCONTROL) + 8140 FORMAT(' x'' mom. Cl''|',20(A,' Cld',I2.2,' =',F11.6)) C WRITE(LU,8150) (' ',K, CMTOT_D(K), K=1, NCONTROL) - 8150 FORMAT(' y mom. Cm |' ,20(A,' Cmd',I1,' =',F11.6)) + 8150 FORMAT(' y mom. Cm |' ,20(A,' Cmd',I2.2,' =',F11.6)) C - WRITE(LU,8160) (' ',K,CNTOT_D(K), K=1, NCONTROL) - 8160 FORMAT(' z'' mom. Cn''|',20(A,' Cnd',I1,' =',F11.6)) + WRITE(LU,8160) (' ',K,DIR*CNTOT_D(K), K=1, NCONTROL) + 8160 FORMAT(' z'' mom. Cn''|',20(A,' Cnd',I2.2,' =',F11.6)) C WRITE(LU,8170) (' ',K, CDFF_D(K), K=1, NCONTROL) - 8170 FORMAT(' Trefftz drag|' ,20(A,'CDffd',I1,' =',F11.6)) + 8170 FORMAT(' Trefftz drag|' ,20(A,'CDffd',I2.2,' =',F11.6)) C WRITE(LU,8180) (' ',K, SPANEF_D(K), K=1, NCONTROL) - 8180 FORMAT(' span eff. |' ,20(A,' ed',I1,' =',F11.6)) + 8180 FORMAT(' span eff. |' ,20(A,' ed',I2.2,' =',F11.6)) C WRITE(LU,*) WRITE(LU,*) @@ -1011,30 +1273,30 @@ SUBROUTINE DERMATS(LU) IF(NDESIGN.GT.0) THEN C WRITE(LU,8206) (GNAME(K), K, K=1, NDESIGN) - 8206 FORMAT(/14X,20(4X,A12, ' g',I1,' ')) + 8206 FORMAT(/14X,20(4X,A12, ' g',I2.2,' ')) WRITE(LU,8207) (' ',K=1, NDESIGN) 8207 FORMAT( 14X,20(3X,A,'----------------')) C WRITE(LU,8210) (' ',K,CLTOT_G(K), K=1, NDESIGN) - 8210 FORMAT(' z'' force CL |' ,20(A,' CLg',I1,' =',F11.6)) + 8210 FORMAT(' z'' force CL |' ,20(A,' CLg',I2.2,' =',F11.6)) C WRITE(LU,8220) (' ',K,CYTOT_G(K), K=1, NDESIGN) - 8220 FORMAT(' y force CY |' ,20(A,' CYg',I1,' =',F11.6)) + 8220 FORMAT(' y force CY |' ,20(A,' CYg',I2.2,' =',F11.6)) C WRITE(LU,8230) (' ',K,DIR*CRTOT_G(K), K=1, NDESIGN) - 8230 FORMAT(' x'' mom. Cl''|' ,20(A,' Clg',I1,' =',F11.6)) + 8230 FORMAT(' x'' mom. Cl''|' ,20(A,' Clg',I2.2,' =',F11.6)) C WRITE(LU,8240) (' ',K, CMTOT_G(K), K=1, NDESIGN) - 8240 FORMAT(' y mom. Cm |' ,20(A,' Cmg',I1,' =',F11.6)) + 8240 FORMAT(' y mom. Cm |' ,20(A,' Cmg',I2.2,' =',F11.6)) C WRITE(LU,8250) (' ',K,DIR*CNTOT_G(K), K=1, NDESIGN) - 8250 FORMAT(' z'' mom. Cn''|',20(A,' Cng',I1,' =',F11.6)) + 8250 FORMAT(' z'' mom. Cn''|',20(A,' Cng',I2.2,' =',F11.6)) C WRITE(LU,8260) (' ',K, CDFF_G(K), K=1, NDESIGN) - 8260 FORMAT(' Trefftz drag|',20(A,'CDffg',I1,' =',F11.6)) + 8260 FORMAT(' Trefftz drag|',20(A,'CDffg',I2.2,' =',F11.6)) C WRITE(LU,8270) (' ',K, SPANEF_G(K), K=1, NDESIGN) - 8270 FORMAT(' span eff. |',20(A,' eg',I1,' =',F11.6)) + 8270 FORMAT(' span eff. |',20(A,' eg',I2.2,' =',F11.6)) C WRITE(LU,*) WRITE(LU,*) @@ -1053,21 +1315,19 @@ SUBROUTINE DERMATS(LU) 8402 FORMAT(/' Clb Cnr / Clr Cnb =', F11.6, & ' ( > 1 if spirally stable )') ENDIF - C RETURN END ! DERMATS - - SUBROUTINE DERMATB(LU) + SUBROUTINE DERMATB(LU, USEMRF) C--------------------------------------------------------- C Calculates and outputs stability derivative matrix C in body axes C--------------------------------------------------------- INCLUDE 'AVL.INC' + LOGICAL USEMRF CHARACTER*50 SATYPE - REAL WROT_RX(3), WROT_RZ(3), WROT_A(3) C CALL GETSA(LNASA_SA,SATYPE,DIR) C @@ -1077,6 +1337,131 @@ SUBROUTINE DERMATB(LU) C---- calculate forces and sensitivities CALL AERO C +C... machine-readable output format + IF(USEMRF) THEN + CALL MRFTOT(LU, 'DERMATB') +C + WRITE(LU,'(/A)') 'Geometry-axis derivatives...' +C + WRITE(LU,'(A)') + & 'axial vel. u, sideslip vel. v, normal vel. w' + WRITE(LU,'(3(ES23.15),2X,A)') + & - CXTOT_U(1), + & -DIR*CXTOT_U(2), + & - CXTOT_U(3), + & '| x force CX : CXu, CXv, CXw' + WRITE(LU,'(3(ES23.15),2X,A)') + & -DIR*CYTOT_U(1), + & - CYTOT_U(2), + & -DIR*CYTOT_U(3), + & '| y force CY : CYu, CYv, CYw' + WRITE(LU,'(3(ES23.15),2X,A)') + & - CZTOT_U(1), + & -DIR*CZTOT_U(2), + & - CZTOT_U(3), + & '| z force CZ : CZu, CZv, CZw' + WRITE(LU,'(3(ES23.15),2X,A)') + & - CRTOT_U(1), + & -DIR*CRTOT_U(2), + & - CRTOT_U(3), + & '| x mom. Cl : Clu, Clv, Clw' + WRITE(LU,'(3(ES23.15),2X,A)') + & -DIR*CMTOT_U(1), + & - CMTOT_U(2), + & -DIR*CMTOT_U(3), + & '| y mom. Cm : Cmu, Cmv, Cmw' + WRITE(LU,'(3(ES23.15),2X,A)') + & - CNTOT_U(1), + & -DIR*CNTOT_U(2), + & - CNTOT_U(3), + & '| z mom. Cn : Cnu, Cnv, Cnw' +C + WRITE(LU,'(A)') 'roll rate p, pitch rate q, yaw rate r' + WRITE(LU,'(3(ES23.15),2X,A)') + & CXTOT_U(4)*2.0/BREF, + & DIR*CXTOT_U(5)*2.0/CREF, + & CXTOT_U(6)*2.0/BREF, + & '| x force CX : CXp, CXq ,CXr' + WRITE(LU,'(3(ES23.15),2X,A)') + & DIR*CYTOT_U(4)*2.0/BREF, + & CYTOT_U(5)*2.0/CREF, + & DIR*CYTOT_U(6)*2.0/BREF, + & '| y force CY : CYp, CYq, CYr' + WRITE(LU,'(3(ES23.15),2X,A)') + & CZTOT_U(4)*2.0/BREF, + & DIR*CZTOT_U(5)*2.0/CREF, + & CZTOT_U(6)*2.0/BREF, + & '| z force CZ : CZp, CZq, CZr' + WRITE(LU,'(3(ES23.15),2X,A)') + & CRTOT_U(4)*2.0/BREF, + & DIR*CRTOT_U(5)*2.0/CREF, + & CRTOT_U(6)*2.0/BREF, + & '| x mom. Cl : Clp, Clq, Clr' + WRITE(LU,'(3(ES23.15),2X,A)') + & DIR*CMTOT_U(4)*2.0/BREF, + & CMTOT_U(5)*2.0/CREF, + & DIR*CMTOT_U(6)*2.0/BREF, + & '| y mom. Cm : Cmp, Cmq, Cmr' + WRITE(LU,'(3(ES23.15),2X,A)') + & CNTOT_U(4)*2.0/BREF, + & DIR*CNTOT_U(5)*2.0/CREF, + & CNTOT_U(6)*2.0/BREF, + & '| z mom. Cn : Cnp, Cnq, Cnr' +C + WRITE(LU,'(I4,2X,A)') NCONTROL, '| # control vars' + IF(NCONTROL.GT.0) THEN + DO K=1, NCONTROL + WRITE(LU,'(A)') DNAME(K) + ENDDO + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & (DIR*CXTOT_D(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| x force CX : CXd*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & ( CYTOT_D(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| y force CY : CYd*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & (DIR*CZTOT_D(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| z force CZ : CZd*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & (DIR*CRTOT_D(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| x mom. Cl : Cld*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & ( CMTOT_D(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| y mom. Cm : Cmd*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & (DIR*CNTOT_D(K), K=1, NCONTROL) + WRITE(LU,'(2X,A)') '| z mom. Cn : Cnd*' + ENDIF +C + WRITE(LU,'(I4,2X,A)') NDESIGN, '| # design vars' + IF(NDESIGN.GT.0) THEN + DO K=1, NDESIGN + WRITE(LU,'(A)') GNAME(K) + ENDDO + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & (DIR*CXTOT_G(K), K=1, NDESIGN) + WRITE(LU,'(2X,A)') '| x force CX : CXg*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & ( CYTOT_G(K), K=1, NDESIGN) + WRITE(LU,'(2X,A)') '| y force CY : CYg*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & (DIR*CZTOT_G(K), K=1, NDESIGN) + WRITE(LU,'(2X,A)') '| z force CZ : CZg*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & (DIR*CRTOT_G(K), K=1, NDESIGN) + WRITE(LU,'(2X,A)') '| x mom. Cl : Clg*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & ( CMTOT_G(K), K=1, NDESIGN) + WRITE(LU,'(2X,A)') '| y mom. Cm : Cmg*' + WRITE(LU,'(40(ES23.15))',ADVANCE='NO') + & (DIR*CNTOT_G(K), K=1, NDESIGN) + WRITE(LU,'(2X,A)') '| z mom. Cn : Cng*' + ENDIF +C + RETURN + ENDIF +C +C... original output CALL OUTTOT(LU) C WRITE(LU,7004) @@ -1187,27 +1572,27 @@ SUBROUTINE DERMATB(LU) IF(NCONTROL.GT.0) THEN C WRITE(LU,8106) (DNAME(K), K, K=1, NCONTROL) - 8106 FORMAT(/14X,20(4X,A12, ' d',I1,' ')) + 8106 FORMAT(/14X,20(4X,A12, ' d',I2.2,' ')) WRITE(LU,8107) (' ',K=1, NCONTROL) 8107 FORMAT( 14X,20(3X,A,'----------------')) C WRITE(LU,8110) (' ',K,DIR*CXTOT_D(K), K=1, NCONTROL) - 8110 FORMAT(' x force CX |',20(A,' CXd',I1,' =',F11.6)) + 8110 FORMAT(' x force CX |',20(A,' CXd',I2.2,' =',F11.6)) C WRITE(LU,8120) (' ',K, CYTOT_D(K), K=1, NCONTROL) - 8120 FORMAT(' y force CY |',20(A,' CYd',I1,' =',F11.6)) + 8120 FORMAT(' y force CY |',20(A,' CYd',I2.2,' =',F11.6)) C WRITE(LU,8130) (' ',K,DIR*CZTOT_D(K), K=1, NCONTROL) - 8130 FORMAT(' z force CZ |',20(A,' CZd',I1,' =',F11.6)) + 8130 FORMAT(' z force CZ |',20(A,' CZd',I2.2,' =',F11.6)) C WRITE(LU,8140) (' ',K,DIR*CRTOT_D(K), K=1, NCONTROL) - 8140 FORMAT(' x mom. Cl |',20(A,' Cld',I1,' =',F11.6)) + 8140 FORMAT(' x mom. Cl |',20(A,' Cld',I2.2,' =',F11.6)) C WRITE(LU,8150) (' ',K, CMTOT_D(K), K=1, NCONTROL) - 8150 FORMAT(' y mom. Cm |',20(A,' Cmd',I1,' =',F11.6)) + 8150 FORMAT(' y mom. Cm |',20(A,' Cmd',I2.2,' =',F11.6)) C WRITE(LU,8160) (' ',K,DIR*CNTOT_D(K), K=1, NCONTROL) - 8160 FORMAT(' z mom. Cn |',20(A,' Cnd',I1,' =',F11.6)) + 8160 FORMAT(' z mom. Cn |',20(A,' Cnd',I2.2,' =',F11.6)) C WRITE(LU,*) WRITE(LU,*) @@ -1217,27 +1602,27 @@ SUBROUTINE DERMATB(LU) IF(NDESIGN.GT.0) THEN C WRITE(LU,8206) (GNAME(K), K, K=1, NDESIGN) - 8206 FORMAT(/14X,20(4X,A12, ' g',I1,' ')) + 8206 FORMAT(/14X,20(4X,A12, ' g',I2.2,' ')) WRITE(LU,8207) (' ',K=1, NDESIGN) 8207 FORMAT( 14X,20(3X,A,'----------------')) C WRITE(LU,8210) (' ',K,DIR*CXTOT_G(K), K=1, NDESIGN) - 8210 FORMAT(' x force CX |',20(A,' CXg',I1,' =',F11.6)) + 8210 FORMAT(' x force CX |',20(A,' CXg',I2.2,' =',F11.6)) C WRITE(LU,8220) (' ',K, CYTOT_G(K), K=1, NDESIGN) - 8220 FORMAT(' y force CY |',20(A,' CYg',I1,' =',F11.6)) + 8220 FORMAT(' y force CY |',20(A,' CYg',I2.2,' =',F11.6)) C WRITE(LU,8230) (' ',K,DIR*CZTOT_G(K), K=1, NDESIGN) - 8230 FORMAT(' z force CZ |',20(A,' CZg',I1,' =',F11.6)) + 8230 FORMAT(' z force CZ |',20(A,' CZg',I2.2,' =',F11.6)) C WRITE(LU,8240) (' ',K,DIR*CRTOT_G(K), K=1, NDESIGN) - 8240 FORMAT(' x mom. Cl |',20(A,' Clg',I1,' =',F11.6)) + 8240 FORMAT(' x mom. Cl |',20(A,' Clg',I2.2,' =',F11.6)) C WRITE(LU,8250) (' ',K, CMTOT_G(K), K=1, NDESIGN) - 8250 FORMAT(' y mom. Cm |',20(A,' Cmg',I1,' =',F11.6)) + 8250 FORMAT(' y mom. Cm |',20(A,' Cmg',I2.2,' =',F11.6)) C WRITE(LU,8260) (' ',K,DIR*CNTOT_G(K), K=1, NDESIGN) - 8260 FORMAT(' z mom. Cn |',20(A,' Cng',I1,' =',F11.6)) + 8260 FORMAT(' z mom. Cn |',20(A,' Cng',I2.2,' =',F11.6)) C WRITE(LU,*) WRITE(LU,*) diff --git a/src/atime.f b/src/atime.f new file mode 100644 index 0000000..132e247 --- /dev/null +++ b/src/atime.f @@ -0,0 +1,656 @@ +C*********************************************************************** +C Module: atime.f +C +C Copyright (C) 2002 Mark Drela, Harold Youngren +C +C This program is free software; you can redistribute it and/or modify +C it under the terms of the GNU General Public License as published by +C the Free Software Foundation; either version 2 of the License, or +C (at your option) any later version. +C +C This program is distributed in the hope that it will be useful, +C but WITHOUT ANY WARRANTY; without even the implied warranty of +C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +C GNU General Public License for more details. +C +C You should have received a copy of the GNU General Public License +C along with this program; if not, write to the Free Software +C Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. +C*********************************************************************** +C +C NOT USED IN CURRENT AVL CODE!!! +C + SUBROUTINE TIME +C--------------------------------------------- +C Time-domain integration driver +C--------------------------------------------- + INCLUDE 'AVL.INC' + INCLUDE 'AVLPLT.INC' + LOGICAL ERROR, LOK, LWRIT, SAVED, OVERLAY +C + CHARACTER*1 ITEM, ANS, CHKEY + CHARACTER*2 OPT, CHPLOT + CHARACTER*4 COMAND, ITEMC + CHARACTER*256 FNOUT, FNNEW, FNSYS, FNVB + CHARACTER*80 LINE, COMARG, PROMPT, RTNEW + CHARACTER SATYPE*50, ROTTYPE*50 +C + REAL*8 ASYS(JETOT,JETOT),BSYS(JETOT,NDMAX),RSYS(JETOT) +C + PARAMETER (JTMAX=JETOT+NDMAX) + REAL AMAT(JTMAX,JTMAX),RMAT(JTMAX),WORK(JTMAX) + REAL DTPAR(JEMAX), DTCON(NDMAX) + INTEGER IMAT(JTMAX) +C + REAL RINPUT(20), RINP(20) + INTEGER IINPUT(20), IINP(20) +C + IRUN0 = IRUN +C + FNVB = ' ' +C + CALL GETSA(LNASA_SA,SATYPE,DIR) +C + IF(LSA_RATES) THEN + ROTTYPE = 'Rates,moments about Stability axes' + ELSE + ROTTYPE = 'Rates,moments about Body axes' + ENDIF +C + 1000 FORMAT (A) +C +C---- force computation of plot limits + TMIN = 0. + TMAX = 0. + FMIN = 0. + FMAX = 0. +C +C + LPLOT = .FALSE. + LWRIT = .FALSE. +C + FNOUT = ' ' +C + CHPLOT = ' ' +C +c EYEPTZ = 180. +c EYEPTY = 90. +c EYEPTX = 0. +c ROBINV = 0. +C + OVERLAY = .FALSE. +C +C================================================================= +C---- start of user interaction loop + 800 CONTINUE +C + WRITE(*,1050) + 1050 FORMAT( + & /' Run-case parameters for candidate initial states... ') +C + LU = 6 + CALL RUNLST(LU,IRUN) +C +C +C + WRITE(*,1052) + 1052 FORMAT( + & ' ==========================================================' + & //' "#" select run case for initial state' + & /' M odify parameters' + & //' eX ecute time march calculation' + & //' V iew time march calculation' + & /' P lot time traces' + & /' W rite time traces to screen or file' + & /' D ata file overlay toggle' + & //' Z oom' + & /' U nzoom') +C +C A B C D E F G H I J K L M N O P Q R S T U V W X Y Z +C x x x x x x x x + + 810 CONTINUE + CALL ASKC(' .TIME^',COMAND,COMARG) +C + 815 CONTINUE +C------------------------------------------------------ + IF (COMAND.EQ.' ') THEN + IF(LPLOT) CALL PLEND + LPLOT = .FALSE. + CALL CLRZOOM + IRUN = IRUN0 + RETURN +C + ELSEIF(COMAND.EQ.'? ') THEN + GO TO 800 +C + ENDIF +C +C------------------------------------------------------------------- +C---- see if command is an integer + NINPUT = 1 + CALL GETINT(COMAND,IINPUT,NINPUT,ERROR) + IF(.NOT.ERROR .AND. NINPUT.GE.1 + & .AND. COMAND(1:1).NE.'T' + & .AND. COMAND(1:1).NE.'F' ) THEN +C----- command is an integer... new case index? + IF(IINPUT(1).LT.0 .OR. IINPUT(1).GT.NRUN) THEN + WRITE(*,*) + WRITE(*,*) '* Selected new run case is not defined' + GO TO 800 + ELSE +C------ valid new run case selected... go back to menu + IRUN = IINPUT(1) + GO TO 800 + ENDIF + ENDIF +C +C------------------------------------------------------------------- +C---- extract command line numeric arguments + DO I=1, 20 + IINPUT(I) = 0 + RINPUT(I) = 0.0 + ENDDO + NINPUT = 20 + CALL GETINT(COMARG,IINPUT,NINPUT,ERROR) + NINPUT = 20 + CALL GETFLT(COMARG,RINPUT,NINPUT,ERROR) +C +C----------------------------------------------------------------------- + IF (COMAND.EQ.'M ') THEN + CALL PARMOD(IRUN) +C +C------------------------------------------------------------------- + ELSEIF(COMAND .EQ. 'X ') THEN +C------ execute time march calculation +C + 22 RINPUT(1) = DELTAT + RINPUT(2) = FLOAT(NTSTEPS) + WRITE(*,2022) DELTAT, NTSTEPS + 2022 FORMAT( + & /' Enter delta(t), Nstep (-Nstep to append steps)', F9.5, I8) + CALL READR(2,RINPUT,ERROR) + IF(ERROR) GO TO 22 +C +C------ abort? + IF(RINPUT(1).LE.0.0 .OR. RINPUT(2).EQ.0.0) GO TO 510 +C + DELTAT = RINPUT(1) + NTSTEPS = INT( RINPUT(2) + SIGN(0.01,RINPUT(2)) ) + NTSTEP = IABS(NTSTEPS) +C +C------ append to current time trace? + LTAPP = NTSTEPS .LT. 0 +C +C------ set up marching indices and initial state + IF(LTAPP) THEN +C------- make sure we have a valid time trace up to now + IF(NTLEV.LT.2) THEN +C--------- bug out + WRITE(*,*) 'No time march available to append to.' + GO TO 510 + ENDIF +C +C------- will just add NTSTEP points starting at current time level ITLEV +C- (any levels after ITLEV will be overwritten) + ITLEV1 = ITLEV + ITLEV2 = ITLEV + NTSTEP +C + ELSE +C------- will start new sequence from current run case IRUN + ITLEV1 = 1 + ITLEV2 = 1 + NTSTEP +C + IR = IRUN +C + CALL RUNCHK(IR,LOK) + IF(.NOT.LOK) THEN + WRITE(*,*) '** No initial state. Ill-posed run case', IR + GO TO 100 + ENDIF +C +C------- converge initial state + NITER = 10 + INFO = 0 + CALL EXEC(NITER,INFO,IR) +C + VEE = PARVAL(IPVEE,IR) + ROT = VEE/UNITL +C + RHO = PARVAL(IPRHO,IR) + QS = 0.5*RHO*VEE**2 * SREF*UNITL**2 +C + WEIGHT = PARVAL(IPMASS,IR)*PARVAL(IPGEE,IR) +C +C------- set initial state + IT = ITLEV1 + TLEV(IT) = 0. + TPARS(KPBANK,IT) = PARVAL(IPPHI,IR)*DTR + TPARS(KPELEV,IT) = PARVAL(IPTHE,IR)*DTR + TPARS(KPHEAD,IT) = PARVAL(IPPSI,IR)*DTR + TPARS(KPVINF,IT) = PARVAL(IPVEE ,IR) + TPARS(KPALFA,IT) = PARVAL(IPALFA,IR)*DTR + TPARS(KPBETA,IT) = PARVAL(IPBETA,IR)*DTR + TPARS(KPCLIF,IT) = CLTOT + TPARS(KPCDRG,IT) = CDTOT + TPARS(KPCDIN,IT) = CDITOT + TPARS(KPCMOM,IT) = CMTOT + TPARS(KPLIFT,IT) = CLTOT*QS + TPARS(KPDRAG,IT) = CDTOT*QS + TPARS(KPMACH,IT) = 0.0 + TPARS(KPALTK,IT) = 0.0 + TPARS(KPDENS,IT) = RHO + TPARS(KPVSOU,IT) = 1.0 + TPARS(KPLOAD,IT) = CLTOT*QS / WEIGHT +C + TPARV(1,KPVEL,IT) = VINF(1)*VEE + TPARV(2,KPVEL,IT) = -VINF(2)*VEE + TPARV(3,KPVEL,IT) = VINF(3)*VEE + TPARV(1,KPROT,IT) = -WROT(1)*ROT + TPARV(2,KPROT,IT) = WROT(2)*ROT + TPARV(3,KPROT,IT) = -WROT(3)*ROT + TPARV(1,KPPOS,IT) = 0.0 + TPARV(2,KPPOS,IT) = 0.0 + TPARV(3,KPPOS,IT) = 0.0 + TPARV(1,KPFOR,IT) = CXTOT + TPARV(2,KPFOR,IT) = CYTOT + TPARV(3,KPFOR,IT) = CZTOT + TPARV(1,KPMOM,IT) = CRTOT + TPARV(2,KPMOM,IT) = CMTOT + TPARV(3,KPMOM,IT) = CNTOT +C + DO N = 1, NCONTROL + TPARD(N,IT) = DELCON(N) + ENDDO + ENDIF +C +C + DO 50 ITLEV = ITLEV1+1, ITLEV2 + IT = ITLEV + ITM = ITLEV - 1 + ITL = MAX( 1 , ITLEV - 2 ) +C +C-------- initial guess for current time from previous time + CALL TCOPY(ITLEV-1,ITLEV) +C + TLEV(ITLEV) = TLEV(ITLEV-1) + DELTAT +C +cC-------- better initial guess by extrapolating in time +c IF(IPTIME.GT.2) CALL QPRED(IPTIME) +c PSTEADY(IPTIME) = .FALSE. +C +C-------- set time-differencing weights + CALL TDSET(ITLEV,TLEV,TDER) +C + WRITE(*,*) + WRITE(*,*) + & 'Converging time level',ITLEV,' ... t = ',TLEV(ITLEV) +C + + + INFO = 1 + ETOL = 1.0E-5 + IF(COMAND .EQ. 'X ') THEN + CALL SYSMAT(IR,ASYS,BSYS,RSYS,NSYS) + ELSE + CALL APPMAT(IR,ASYS,BSYS,RSYS,NSYS) + ENDIF +C +cc LU = 6 +cc CALL SYSSHO(LU,ASYS,BSYS,RSYS,NSYS) +C + NMAT = JETOT + NCONTROL + DO K = 1, NMAT + RMAT(K) = 0. + DO J = 1, NMAT + AMAT(K,J) = 0. + ENDDO + ENDDO +C +C-------- set total residuals and Jacobian + DO K = 1, JETOT + RMAT(K) = RSYS(K) + & - TDER(1)*TPAR(K,IT ) + & - TDER(2)*TPAR(K,ITM) + & - TDER(3)*TPAR(K,ITL) + DO J = 1, JETOT + AMAT(K,J) = ASYS(K,J) - TDER(1) + ENDDO + DO L = 1, NCONTROL + J = JETOT + L + AMAT(K,J) = BSYS(K,L) + ENDDO + ENDDO +C +C-------- set control-variable residuals and Jacobian + DO N = 1, NCONTROL + K = JETOT + N + RMAT(K) = TCON(N,IT) - TCONSP(N,IT) + AMAT(K,K) = 1.0 + ENDDO +C + CALL LUDCMP(JTMAX,NMAT,AMAT,IMAT,WORK) + CALL BAKSUB(JTMAX,NMAT,AMAT,IMAT,RMAT) +C + DRMS = 0. + DMAX = 0. + JMAX = 0 + DO J = 1, JETOT + DTPAR(J) = -RMAT(J) +C + IF(ABS(DTPAR(J))/TREF(J) .GT. ABS(DMAX)) THEN + DMAX = DTPAR(J) + JMAX = J + ENDIF + DRMS = DRMS + (DTPAR(J)/TREF(J))**2 + ENDDO + DO N = 1, NCONTROL + J = JETOT + N + DTCON(N) = -RMAT(J) +C + IF(ABS(DTCON(N))/TREF(J) .GT. ABS(DMAX)) THEN + DMAX = DTCON(N) + JMAX = J + ENDIF + DRMS = DRMS + (DTCON(N)/TREF(J))**2 + ENDDO +C + DRMS = SQRT( DRMS/FLOAT(JETOT+NCONTROL) ) +C +C + RLX = 1.0 + DO J = 1, JETOT + TPAR(J,IT) = TPAR(J,IT) + RLX*DTPAR(J) + ENDDO + DO N = 1, NCONTROL + TCON(N,IT) = TCON(N,IT) + RLX*DTCON(N) + ENDDO +C + VSQ = TPAR(JUE,IT)**2 + TPAR(JVE,IT)**2 + TPAR(JWE,IT)**2 + VEE = SQRT(VSQ) + ROT = VEE/UNITL +C + IF(VEE .EQ. 0.0) THEN + VINV = 1.0 + RINV = 1.0 + ELSE + VINV = 1.0/VEE + RINV = 1.0/ROT + ENDIF +C + VINF(1) = TPAR(JEU ,IT)*VINV + VINF(2) = -TPAR(JEV ,IT)*VINV + VINF(3) = TPAR(JEW ,IT)*VINV + WROT(1) = -TPAR(JEP ,IT)*RINV + WROT(2) = TPAR(JEQ ,IT)*RINV + WROT(3) = -TPAR(JER ,IT)*RINV +C + V13 = SQRT( VINF(1)**2 + VINF(3)**2 ) +C + ALFA = ATAN2( VINF(3) , VINF(1) ) + BETA = ATAN2( -VINF(2) , V13 ) +C + PARVAL(IPALFA,IR) = ALFA/DTR + PARVAL(IPBETA,IR) = BETA/DTR + PARVAL(IPROTX,IR) = WROT(1)*0.5*BREF + PARVAL(IPROTY,IR) = WROT(2)*0.5*CREF + PARVAL(IPROTZ,IR) = WROT(3)*0.5*BREF + PARVAL(IPCL ,IR) = CLTOT +C + PARVAL(IPALFA,IR) = ALFA + PARVAL(IPBETA,IR) = BETA + PARVAL(IPROTX,IR) = WROT(1) + PARVAL(IPROTY,IR) = WROT(2) + PARVAL(IPROTZ,IR) = WROT(3) + PARVAL(IPPHI,IR) = TPAR(JEPH,IT) + PARVAL(IPTHE,IR) = TPAR(JETH,IT) + PARVAL(IPPSI,IR) = TPAR(JEPS,IT) +C +C-------- sum AIC matrices to get GAM,SRC,DBL + CALL GAMSUM +C +C-------- sum AIC matrices to get WC,WV + CALL VELSUM +C +C-------- compute forces + CALL AERO +C +C + IF(DRMS .LT. EPS .AND. ABS(DMAX) .LT. 10.0*EPS) THEN + GO TO 100 + ENDIF + + +C + 100 CONTINUE +C + CHPLOT = 'P ' +C +C------------------------------------------------------------------- +C---- Blowup window + ELSEIF(COMAND.EQ.'B ') THEN + IF(INDEX('P',CHPLOT(1:1)) .EQ. 0) THEN + WRITE(*,*) 'No plot to blow up' + GO TO 810 + ENDIF +C + DCROSS = 2.0 + WRITE(*,*) 'Mark off corners of blowup area' + CALL GETCURSORXY(XC1,YC1,CHKEY) + CALL PLSYMB(XC1,YC1,DCROSS,3,0.0,0) + CALL PLFLUSH + CALL GETCURSORXY(XC2,YC2,CHKEY) + CALL PLSYMB(XC2,YC2,DCROSS,3,0.0,0) + CALL PLFLUSH +C + XE1 = TMIN+(XC1-XORG0)/TFAC + YE1 = FMIN+(YC1-YORG0)/FFAC + XE2 = TMIN+(XC2-XORG0)/TFAC + YE2 = FMIN+(YC2-YORG0)/FFAC +C + TMIN = MIN( XE1 , XE2 ) + TMAX = MAX( XE1 , XE2 ) + FMIN = MIN( YE1 , YE2 ) + FMAX = MAX( YE1 , YE2 ) +C + CALL AXISADJ(TMIN,TMAX, TTOT, TDEL, NANN) + CALL AXISADJ(FMIN,FMAX, FTOT, FDEL, NANN) +C + COMAND = 'P ' + GO TO 815 +C +C------------------------------------------------------------------- +C---- Normal size + ELSEIF(COMAND.EQ.'R ') THEN + TMIN = 0. + TMAX = 0. + FMIN = 0. + FMAX = 0. +C + COMAND = 'P ' + GO TO 815 +C +C------------------------------------------------------------------- +C---- Annotate + ELSEIF(COMAND.EQ.'A ') THEN + IF(LPLOT) THEN + CALL ANNOT(CH) + ELSE + WRITE(*,*) 'No active plot' + ENDIF +C +C------------------------------------------------------------------- +C---- Hardcopy + ELSEIF(COMAND.EQ.'H ') THEN + IF(LPLOT) CALL PLEND + LPLOT = .FALSE. + CALL REPLOT(IDEVH) +C +C------------------------------------------------------------------- +C---- write eigenvalues + ELSEIF(COMAND.EQ.'W ') THEN + IF(FEVDEF(1:1).EQ.' ') THEN +C------ set default filename + KDOT = INDEX(FILDEF,'.') + IF(KDOT.EQ.0) THEN + CALL SLEN(FILDEF,NFIL) + FEVDEF = FILDEF(1:NFIL) // '.eig' + ELSE + FEVDEF = FILDEF(1:KDOT) // 'eig' + ENDIF + ENDIF +C + CALL SLEN(FEVDEF,NFE) + NFE = MAX( NFE , 1 ) + WRITE(*,2040) FEVDEF(1:NFE) + 2040 FORMAT(' Enter eigenvalue save filename: ', A) + READ (*,1000) FNNEW + IF(FNNEW.NE.' ') FEVDEF = FNNEW + CALL EIGOUT(FEVDEF, IRUN1,IRUN2, SAVED) +C +C------------------------------------------------------------------- +C---- toggle overlay flag + ELSEIF(COMAND.EQ.'D ') THEN + OVERLAY = .NOT.OVERLAY +C + IF(OVERLAY) THEN + IF(FEVDEF(1:1).EQ.' ') THEN +C------- set default filename + KDOT = INDEX(FILDEF,'.') + IF(KDOT.EQ.0) THEN + CALL SLEN(FILDEF,NFIL) + FEVDEF = FILDEF(1:NFIL) // '.eig' + ELSE + FEVDEF = FILDEF(1:KDOT) // 'eig' + ENDIF + ENDIF +C + CALL SLEN(FEVDEF,NFE) + NFE = MAX( NFE , 1 ) + WRITE(*,2050) FEVDEF(1:NFE) + 2050 FORMAT(' Enter eigenvalue data filename: ', A) + READ (*,1000) FNNEW + IF(FNNEW.NE.' ') FEVDEF = FNNEW +C +C------ read the data file + CALL EIGINP(FEVDEF, ERROR) + IF(ERROR) THEN + OVERLAY = .FALSE. + GO TO 800 + ENDIF +C + OVERLAY = .TRUE. + WRITE(*,*) 'Overlay plotting enabled' +C + IF(INDEX('P',CHPLOT(1:1)).NE.0) THEN + COMAND = CHPLOT + GO TO 815 + ENDIF +C + ELSE + WRITE(*,*) 'Overlay plotting disabled' +C + ENDIF +C +C------------------------------------------------------------------- +C---- Zoom in on plot + ELSEIF(COMAND.EQ.'Z ') THEN + IF(INDEX('P',CHPLOT(1:1)) .EQ. 0) THEN + WRITE(*,*) 'No plot to zoom' + GO TO 810 + ENDIF +C + CALL USETZOOM(.TRUE.,.TRUE.) + CALL REPLOT(IDEV) +C +C------------------------------------------------------------------- +C---- Reset zoom on plot + ELSEIF(COMAND.EQ.'U ') THEN + CALL CLRZOOM + CALL REPLOT(IDEV) +C +C------------------------------------------------------------------- + ELSE + WRITE(*,*) + WRITE(*,*) '* Option not recognized' +C + ENDIF + GO TO 800 +C + END ! TIME + + + + SUBROUTINE TCOPY(IT1,IT2) + INCLUDE 'AVL.INC' +C + TLEV(IT2) = TLEV(IT1) + DO J = 1, JETOT + TPAR(J,IT2) = TPAR(J,IT1) + ENDDO + DO N = 1, NCONTROL + TCON(N,IT2) = TCON(N,IT1) + ENDDO +C + RETURN + END ! TCOPY + + + + SUBROUTINE TDSET(ITIME,TIME,QTWT) + REAL TIME(*), QTWT(3) +C + REAL QTWT1(3), QTWT2(3) +C +C---- set differencing blending fraction +C +C- F = 0 1st-order backward difference +C- F = 1 2nd-order backward difference +C + DATA F / 1.0 / +ccc DATA F / 0.0 / +C + IF (ITIME.LT.2) THEN + QTWT1(1) = 0. + QTWT1(2) = 0. + QTWT1(3) = 0. + QTWT2(1) = 0. + QTWT2(2) = 0. + QTWT2(3) = 0. +C + ELSEIF(ITIME.EQ.2) THEN + T0 = TIME(ITIME ) + T1 = TIME(ITIME-1) +C + QTWT1(1) = 1.0/(T0-T1) + QTWT1(2) = -1.0/(T0-T1) + QTWT1(3) = 0. +C + QTWT2(1) = 1.0/(T0-T1) + QTWT2(2) = -1.0/(T0-T1) + QTWT2(3) = 0. +C + ELSE + T0 = TIME(ITIME ) + T1 = TIME(ITIME-1) + T2 = TIME(ITIME-2) +C + QTWT1(1) = 1.0/(T0-T1) + QTWT1(2) = -1.0/(T0-T1) + QTWT1(3) = 0. +C + QTWT2(1) = 1.0/(T0-T1) + 1.0/(T0-T2) + QTWT2(2) = -1.0/(T0-T1) - 1.0/(T0-T2) + & - (T0-T1)/((T0-T2)*(T1-T2)) + QTWT2(3) = (T0-T1)/((T0-T2)*(T1-T2)) +C + ENDIF +C +C---- blend 1st and 2nd order differences + QTWT(1) = (1.0-F)*QTWT1(1) + F*QTWT2(1) + QTWT(2) = (1.0-F)*QTWT1(2) + F*QTWT2(2) + QTWT(3) = (1.0-F)*QTWT1(3) + F*QTWT2(3) +C + RETURN + END ! TDSET + diff --git a/src/atpforc.f b/src/atpforc.f index ab58b40..4a41e52 100644 --- a/src/atpforc.f +++ b/src/atpforc.f @@ -281,33 +281,38 @@ SUBROUTINE TPFORC C 30 CONTINUE C -C...Trefftz-plane drag is kinetic energy in crossflow DWWAKE(JC) = -(NY*VY + NZ*VZ) C - CLFF = CLFF + 2.0*GAMS(JC)* DYT /SREF - CYFF = CYFF - 2.0*GAMS(JC)* DZT /SREF - CDFF = CDFF + GAMS(JC)*(DZT*VY - DYT*VZ)/SREF - DO N = 1, NUMAX - CLFF_U(N) = CLFF_U(N) + 2.0*GAMS_U(JC,N)*DYT/SREF - CYFF_U(N) = CYFF_U(N) - 2.0*GAMS_U(JC,N)*DZT/SREF - CDFF_U(N) = CDFF_U(N) +C...Trefftz-plane drag is kinetic energy in crossflow +C + ISURF = NSURFS(JC) + IF(LFLOAD(ISURF)) THEN +C-------add load from this strip only if it contributes to total load + CLFF = CLFF + 2.0*GAMS(JC)* DYT /SREF + CYFF = CYFF - 2.0*GAMS(JC)* DZT /SREF + CDFF = CDFF + GAMS(JC)*(DZT*VY - DYT*VZ)/SREF + DO N = 1, NUMAX + CLFF_U(N) = CLFF_U(N) + 2.0*GAMS_U(JC,N)*DYT/SREF + CYFF_U(N) = CYFF_U(N) - 2.0*GAMS_U(JC,N)*DZT/SREF + CDFF_U(N) = CDFF_U(N) & + ( GAMS_U(JC,N)*(DZT*VY - DYT*VZ ) & + GAMS(JC) *(DZT*VY_U(N) - DYT*VZ_U(N)) )/SREF - ENDDO - DO N = 1, NCONTROL - CLFF_D(N) = CLFF_D(N) + 2.0*GAMS_D(JC,N)*DYT/SREF - CYFF_D(N) = CYFF_D(N) - 2.0*GAMS_D(JC,N)*DZT/SREF - CDFF_D(N) = CDFF_D(N) + ENDDO + DO N = 1, NCONTROL + CLFF_D(N) = CLFF_D(N) + 2.0*GAMS_D(JC,N)*DYT/SREF + CYFF_D(N) = CYFF_D(N) - 2.0*GAMS_D(JC,N)*DZT/SREF + CDFF_D(N) = CDFF_D(N) & + ( GAMS_D(JC,N)*(DZT*VY - DYT*VZ ) & + GAMS(JC) *(DZT*VY_D(N) - DYT*VZ_D(N)) )/SREF - ENDDO - DO N = 1, NDESIGN - CLFF_G(N) = CLFF_G(N) + 2.0*GAMS_G(JC,N)*DYT/SREF - CYFF_G(N) = CYFF_G(N) - 2.0*GAMS_G(JC,N)*DZT/SREF - CDFF_G(N) = CDFF_G(N) + ENDDO + DO N = 1, NDESIGN + CLFF_G(N) = CLFF_G(N) + 2.0*GAMS_G(JC,N)*DYT/SREF + CYFF_G(N) = CYFF_G(N) - 2.0*GAMS_G(JC,N)*DZT/SREF + CDFF_G(N) = CDFF_G(N) & + ( GAMS_G(JC,N)*(DZT*VY - DYT*VZ ) & + GAMS(JC) *(DZT*VY_G(N) - DYT*VZ_G(N)) )/SREF - ENDDO + ENDDO + ENDIF 40 CONTINUE C C---- Double the X,Z forces, zero Y force for a Y symmetric case diff --git a/src/atrim.f b/src/atrim.f index 65899c9..410bd02 100644 --- a/src/atrim.f +++ b/src/atrim.f @@ -37,29 +37,16 @@ SUBROUTINE TRMSET(COMAND,COMARG,COM,CARG) IR = 1 C C---- for trim command, first number in arg string should be part of command - - if(LVERBOSE) then - WRITE(*,*) 'TRIMSET --- YES ' - WRITE(*,*)'COMAND(1:1): ' , COMAND(1:1) - WRITE(*,*)'COMAND(1:1): ' , COMAND(1:1) - endif - IF(COMAND(1:1) .NE. 'C') THEN LMATCH = .FALSE. - RETURN ENDIF C CNUM = COMARG(1:2) C - if(LVERBOSE) then - WRITE(*,*)'CNUM: ' , CNUM - ENDIF - +C IF (CNUM.EQ.'1 ') THEN - KTRIM = 1 - cc CALL TRIM1(IR) LMATCH = .TRUE. C @@ -74,7 +61,7 @@ SUBROUTINE TRMSET(COMAND,COMARG,COM,CARG) C ENDIF C -C WRITE(*,*)'KTRIM: ', KTRIM +C C---- start here with new run case 100 CONTINUE C @@ -96,6 +83,7 @@ SUBROUTINE TRMSET(COMAND,COMARG,COM,CARG) 101 CONTINUE ENDIF C +C WRITE(*,*) C C---- tag this run case with trim type (not converged yet) ITRIM(IR) = -KTRIM @@ -141,7 +129,9 @@ SUBROUTINE TRMSET(COMAND,COMARG,COM,CARG) & CL .GT. 0.0 ) THEN VEE = SQRT( 2.0*RMASS*GEE / (RHO*SREFD*CL*COSP) ) PARVAL(IPVEE,JR) = VEE -C WRITE(*,7500) 'velocity', JR + if (lverbose) then + WRITE(*,7500) 'velocity', JR + endif ENDIF C C-------- set CL ? @@ -149,7 +139,9 @@ SUBROUTINE TRMSET(COMAND,COMARG,COM,CARG) & VEE .GT. 0.0 ) THEN CL = 2.0*RMASS*GEE / (RHO*SREFD*VEE**2*COSP) PARVAL(IPCL,JR) = CL -C WRITE(*,7500) 'CL', JR + if (lverbose) then + WRITE(*,7500) 'CL', JR + endif ENDIF C IF(SINP .EQ. 0.0) THEN @@ -158,15 +150,19 @@ SUBROUTINE TRMSET(COMAND,COMARG,COM,CARG) RAD = VEE**2 * COSP / (GEE * SINP) ENDIF PARVAL(IPRAD,JR) = RAD -C WRITE(*,7500) 'turn radius', JR + if (lverbose) then + WRITE(*,7500) 'turn radius', JR + endif C FAC = 1.0/COSP PARVAL(IPFAC,JR) = FAC -C WRITE(*,7500) 'load factor', JR + if (lverbose) then + WRITE(*,7500) 'load factor', JR + endif C THE = 0. PARVAL(IPTHE,JR) = THE -C cc WRITE(*,7500) 'elevation', JR +cc WRITE(*,7500) 'elevation', JR C C-------- set up CL and rotation rate constraints IF(RAD .GT. 0.0) THEN @@ -214,7 +210,9 @@ SUBROUTINE TRMSET(COMAND,COMARG,COM,CARG) & CL .GT. 0.0 ) THEN RAD = RMASS / (0.5*RHO * SREFD * CL) PARVAL(IPRAD,JR) = RAD -C WRITE(*,7500) 'turn radius', JR + if (lverbose) then + WRITE(*,7500) 'turn radius', JR + endif ENDIF C C-------- set CL ? @@ -222,7 +220,9 @@ SUBROUTINE TRMSET(COMAND,COMARG,COM,CARG) & CL .EQ. 0.0 ) THEN CL = RMASS / (0.5*RHO * SREFD * RAD) PARVAL(IPCL,JR) = CL -C WRITE(*,7500) 'CL', JR + if (lverbose) then + WRITE(*,7500) 'CL', JR + endif ENDIF C C-------- set load factor ? @@ -232,7 +232,9 @@ SUBROUTINE TRMSET(COMAND,COMARG,COM,CARG) & GEE .GT. 0.0 ) THEN FAC = 0.5*RHO*VEE**2 * SREFD * CL / (RMASS*GEE) PARVAL(IPFAC,JR) = FAC -C WRITE(*,7500) 'load factor', JR + if (lverbose) then + WRITE(*,7500) 'load factor', JR + endif ENDIF C C-------- set velocity ? @@ -242,12 +244,14 @@ SUBROUTINE TRMSET(COMAND,COMARG,COM,CARG) & GEE .GT. 0.0 ) THEN VEE = SQRT( FAC * RMASS*GEE / (0.5*RHO*SREFD*CL) ) PARVAL(IPVEE,JR) = VEE -C WRITE(*,7500) 'velocity', JR + if (lverbose) then + WRITE(*,7500) 'velocity', JR + endif ENDIF C THE = 0. PARVAL(IPTHE,JR) = THE -C cc WRITE(*,7500) 'elevation', JR +cc WRITE(*,7500) 'elevation', JR C C-------- set up CL and rotation rate constraints IF(RAD .GT. 0.0) THEN @@ -301,49 +305,45 @@ SUBROUTINE TRMSET(COMAND,COMARG,COM,CARG) C----- jump back here just for menu 10 CONTINUE CALL CFRAC(IR,NRUN,PROMPT,NPR) -C WRITE(*,2000) PROMPT(1:NPR),RTITLE(IR) + if (lverbose) then + WRITE(*,2000) PROMPT(1:NPR),RTITLE(IR) + endif C IF (KTRIM.EQ.1) THEN - if(LVERBOSE) then - WRITE(*,2005) '(level or banked horizontal flight)' - WRITE(*,2105) 'B bank angle = ', PHI , 'deg' - WRITE(*,2105) 'C CL = ', CL , ' ' - WRITE(*,2105) 'V velocity = ', VEE , UNCHV(1:NUV) - WRITE(*,2105) 'M mass = ', RMASS, UNCHM(1:NUM) - WRITE(*,2105) 'D air dens. = ', RHO , UNCHD(1:NUD) - WRITE(*,2105) 'G grav.acc. = ', GEE , UNCHA(1:NUA) - WRITE(*,2105) ' turn rad. = ', RAD , UNCHL(1:NUL) - WRITE(*,2105) ' load fac. = ', FAC , ' ' - WRITE(*,2105) 'X X_cg = ', XCG , 'Lunit' - WRITE(*,2105) 'Y Y_cg = ', YCG , 'Lunit' - WRITE(*,2105) 'Z Z_cg = ', ZCG , 'Lunit' - endif + if (lverbose) then + WRITE(*,2005) '(level or banked horizontal flight)' + WRITE(*,2105) 'B bank angle = ', PHI , 'deg' + WRITE(*,2105) 'C CL = ', CL , ' ' + WRITE(*,2105) 'V velocity = ', VEE , UNCHV(1:NUV) + WRITE(*,2105) 'M mass = ', RMASS, UNCHM(1:NUM) + WRITE(*,2105) 'D air dens. = ', RHO , UNCHD(1:NUD) + WRITE(*,2105) 'G grav.acc. = ', GEE , UNCHA(1:NUA) + WRITE(*,2105) ' turn rad. = ', RAD , UNCHL(1:NUL) + WRITE(*,2105) ' load fac. = ', FAC , ' ' + WRITE(*,2105) 'X X_cg = ', XCG , 'Lunit' + WRITE(*,2105) 'Y Y_cg = ', YCG , 'Lunit' + WRITE(*,2105) 'Z Z_cg = ', ZCG , 'Lunit' + endif C C ELSEIF(KTRIM.EQ.2) THEN - if(LVERBOSE) then - WRITE(*,2005) '(steady pitch rate - looping flight)' - WRITE(*,2105) 'C CL = ', CL , ' ' - WRITE(*,2105) 'V velocity = ', VEE , UNCHV(1:NUV) - WRITE(*,2105) 'M mass = ', RMASS, UNCHM(1:NUM) - WRITE(*,2105) 'D air dens. = ', RHO , UNCHD(1:NUD) - WRITE(*,2105) 'G grav.acc. = ', GEE , UNCHA(1:NUA) - WRITE(*,2105) 'R turn rad. = ', RAD , UNCHL(1:NUL) - WRITE(*,2105) 'L load fac. = ', FAC , ' ' - WRITE(*,2105) 'X X_cg = ', XCG , 'Lunit' - WRITE(*,2105) 'Y Y_cg = ', YCG , 'Lunit' - WRITE(*,2105) 'Z Z_cg = ', ZCG , 'Lunit' + if (lverbose) then + WRITE(*,2005) '(steady pitch rate - looping flight)' + WRITE(*,2105) 'C CL = ', CL , ' ' + WRITE(*,2105) 'V velocity = ', VEE , UNCHV(1:NUV) + WRITE(*,2105) 'M mass = ', RMASS, UNCHM(1:NUM) + WRITE(*,2105) 'D air dens. = ', RHO , UNCHD(1:NUD) + WRITE(*,2105) 'G grav.acc. = ', GEE , UNCHA(1:NUA) + WRITE(*,2105) 'R turn rad. = ', RAD , UNCHL(1:NUL) + WRITE(*,2105) 'L load fac. = ', FAC , ' ' + WRITE(*,2105) 'X X_cg = ', XCG , 'Lunit' + WRITE(*,2105) 'Y Y_cg = ', YCG , 'Lunit' + WRITE(*,2105) 'Z Z_cg = ', ZCG , 'Lunit' endif C C ENDIF C C CALL ASKC(' Enter parameter, value (or # - + N )',COM,CARG) C - -C - -C WRITE(*,*) "COM: ", COM -C WRITE(*,*) "CARG: ", CARG - IF(COM.EQ.' ') THEN C------ just a Return entered... go back RETURN @@ -356,7 +356,7 @@ SUBROUTINE TRMSET(COMAND,COMARG,COM,CARG) C----- add new case after current one C IF(NRUN.EQ.NRMAX) THEN -C WRITE(*,*) + WRITE(*,*) WRITE(*,*) '* Run case array limit NRMAX reached' ELSE NRUN = NRUN + 1 @@ -682,7 +682,9 @@ SUBROUTINE TRMSET(COMAND,COMARG,COM,CARG) ENDIF DO JR = IR1, IR2 PARVAL(IPXCG,JR) = XCG -C WRITE(*,7500) 'x_cg', JR + if (lverbose) then + WRITE(*,7500) 'x_cg', JR + endif ENDDO C C------------------------------------- @@ -694,7 +696,9 @@ SUBROUTINE TRMSET(COMAND,COMARG,COM,CARG) ENDIF DO JR = IR1, IR2 PARVAL(IPYCG,JR) = YCG -C WRITE(*,7500) 'y_cg', JR + if (lverbose) then + WRITE(*,7500) 'y_cg', JR + endif ENDDO C C------------------------------------- @@ -706,7 +710,9 @@ SUBROUTINE TRMSET(COMAND,COMARG,COM,CARG) ENDIF DO JR = IR1, IR2 PARVAL(IPZCG,JR) = ZCG -C WRITE(*,7500) 'z_cg', JR + if (lverbose) then + WRITE(*,7500) 'z_cg', JR + endif ENDDO C C------------------------------------- @@ -718,7 +724,9 @@ SUBROUTINE TRMSET(COMAND,COMARG,COM,CARG) ENDIF DO JR = IR1, IR2 PARVAL(IPCD0,JR) = CD0 -C WRITE(*,7500) 'CDo', JR + if (lverbose) then + WRITE(*,7500) 'CDo', JR + endif ENDDO C C------------------------------------------------------ diff --git a/src/avl.f b/src/avl.f index 9c158d8..cd79a9b 100644 --- a/src/avl.f +++ b/src/avl.f @@ -17,44 +17,8 @@ C along with this program; if not, write to the Free Software C Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. C*********************************************************************** - -C SUBROUTINE AVL - -C INCLUDE 'AVL_test.INC' - - -C DTR = PI/180.0 - -C PI = 3.14 -C DTR = PI/180.0 - -C END - -C C FILE: COMMON_ex.F -C SUBROUTINE FOO -C C REAL PI -C INCLUDE 'AVL_test.INC' - -C C INCLUDE 'AVL_test.INC' -C C COMMON/MASS/PI - -C C C -C C PI = 4.0*ATAN(1.0) -C C PI = 3.14 - -C PRINT*, "PI=",PI -C PRINT*, "DTR=",DTR - -C END -C C END OF COMMON_ex.F -C C PROGRAM AVL - - - - - - +C PROGRAM AVL SUBROUTINE AVL C C======================================================================= C C 3-D Vortex Lattice code. @@ -62,20 +26,18 @@ SUBROUTINE AVL C C See file version_notes.txt for most recent changes. C C======================================================================= INCLUDE 'AVL.INC' -C INCLUDE 'AVLPLT.INC' - LOGICAL ERROR - +C INCLUDE 'AVLPLT.INC' + LOGICAL ERROR, LINPFILE +C CHARACTER*4 COMAND CHARACTER*128 COMARG - CHARACTER*120 FNNEW + CHARACTER*256 FNNEW C REAL RINPUT(20) INTEGER IINPUT(20) C - -C COMMON/CASE_R/PI - VERSION = 3.35 + VERSION = 3.40 C 1000 FORMAT(A) C @@ -92,9 +54,10 @@ SUBROUTINE AVL C PI = 4.0*ATAN(1.0) DTR = PI/180.0 - - - +C +C---- Flag for having valid input data from file + LINPFILE = .FALSE. +C C---- logical units LUINP = 4 ! configuration file LURUN = 7 ! run case file @@ -105,9 +68,7 @@ SUBROUTINE AVL LUSYS = 22 ! dynamic system matrix dump file C C---- set basic defaults - CALL DEFINI - CALL MASINI C C---- initialize Xplot, and AVL plot stuff @@ -125,11 +86,16 @@ SUBROUTINE AVL C C----- no valid geometry... skip reading run and mass files C IF(ERROR) GO TO 100 C C +C C----- This is moved to the loadgeo function +CC----- we have an input dataset to process +C LINPFILE = .TRUE. +CC C C----- set up all parameters C CALL PARSET C C----- process geometry to define strip and vortex data LPLTNEW = .TRUE. + ! TODO remove? CALL ENCALC C C----- initialize state @@ -462,17 +428,12 @@ SUBROUTINE AVL END - SUBROUTINE loadGEO(FILDEF) + SUBROUTINE loadGEO(geom_file) INCLUDE 'AVL.INC' - CHARACTER*120 FILDEF + CHARACTER*256 geom_file LOGICAL ERROR -C LOGICAL LAIC -C LOGICAL LSRD -C LOGICAL LVEL -C LOGICAL LSOL -C LOGICAL LSEN -C LOGICAL LPLTNEW -C +C + FILDEF = geom_file CALL INPUT(LUINP,FILDEF,ERROR) IF(ERROR) THEN WRITE(*,*) @@ -502,13 +463,12 @@ SUBROUTINE loadGEO(FILDEF) LVEL = .FALSE. LSOL = .FALSE. LSEN = .FALSE. -C - + END - SUBROUTINE loadMASS(FMSDEF) + SUBROUTINE loadMASS(mass_file) INCLUDE 'AVL.INC' - CHARACTER*120 FMSDEF + CHARACTER*256 mass_file LOGICAL ERROR C LOGICAL LSOL C LOGICAL LSEN @@ -538,6 +498,7 @@ SUBROUTINE loadMASS(FMSDEF) C CALL STRIP(FMSDEF,NMS) + FMSDEF = mass_file CALL MASGET(LUMAS,FMSDEF,ERROR) IF(ERROR) THEN ELSE @@ -546,14 +507,16 @@ SUBROUTINE loadMASS(FMSDEF) C C CALL MASSHO(6) C CALL APPGET -C WRITE(*,*) -C & '- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -' -C CALL APPSHO(6,RHO0) -C WRITE(*,*) 'RHO0: ', RHO0 -C WRITE(*,*) -C WRITE(*,*) -C & 'Use MSET to apply these mass,inertias to run cases' -C ccc CALL MASPUT(1,NRMAX) + if (LVERBOSE) then + WRITE(*,*) + & '- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -' + CALL APPSHO(6,RHO0) +C + WRITE(*,*) + WRITE(*,*) + & 'Use MSET to apply these mass,inertias to run cases' +ccc CALL MASPUT(1,NRMAX) + endif ENDIF IR1 = 0 @@ -566,7 +529,7 @@ SUBROUTINE loadMASS(FMSDEF) IR2 = IR1 ENDIF C - IF(IR1.LT.1 .OR. IR1.GT.NRUN) WRITE(*,*) "ERROR" + IF(IR1.LT.1 .OR. IR1.GT.NRUN) WRITE(*,*) "ERROR in avl.f" C CALL MASPUT(IR1,IR2) C @@ -574,14 +537,7 @@ SUBROUTINE loadMASS(FMSDEF) LSEN = .FALSE. C - END ! LOAFGEO -C -C SUBROUTINE setALPHA(ANGLE) -C REAL ANGLE - - - ! loadGeo - + END C SUBROUTINE PLINIT C C---- Initialize plotting variables @@ -1078,23 +1034,11 @@ SUBROUTINE RUNINI ICON(IVROTY,IR) = ICROTY ICON(IVROTZ,IR) = ICROTZ C - -c ! if (lverbose)then -c ! WRITE(*,*) "=========================" -c ! end if - C------ default constraint values -c ! if (lverbose)then -c ! WRITE(*,*) "ICALFA", ICALFA -c ! end if DO IC = 1, ICTOT CONVAL(IC,IR) = 0. -C WRITE(*,*) "CONVAL(IC,IR)", CONVAL(IC,IR) -C WRITE(*,*) "IC", IC -C WRITE(*,*) "IR", IR ENDDO -C WRITE(*,*) "=========================" - +C C------ default run case titles RTITLE(IR) = ' -unnamed- ' C @@ -1353,4 +1297,3 @@ SUBROUTINE AOCFIL(FNAME,IFILE) RETURN END - END ! INTE diff --git a/src/build/fileList b/src/build/fileList index 64b98cd..018047e 100644 --- a/src/build/fileList +++ b/src/build/fileList @@ -24,6 +24,8 @@ f77Files =\ spline.f\ userio.f\ eispack.f\ + aoutmrf.f\ + aoml.f\ ad_src/forward_ad_src/aero_d.f\ ad_src/forward_ad_src/aic_d.f\ ad_src/forward_ad_src/cdcl_d.f\ diff --git a/src/f2py/libavl.pyf b/src/f2py/libavl.pyf index 7422acd..773df6b 100644 --- a/src/f2py/libavl.pyf +++ b/src/f2py/libavl.pyf @@ -35,12 +35,12 @@ python module libavl ! in subroutine aero_b ! in :libavl:aero_b.f end subroutine aero_b - subroutine loadgeo(fildef) ! in :libavl:avl.f - character*120 :: fildef + subroutine loadgeo(geom_file) ! in :libavl:avl.f + character*256 :: geom_file end subroutine loadgeo - subroutine loadmass(fmsdef) ! in :libavl:avl.f - character*120 :: fmsdef + subroutine loadmass(mass_file) ! in :libavl:avl.f + character*256 :: mass_file end subroutine loadmass diff --git a/src/includes/AVL.INC b/src/includes/AVL.INC index 002c5d6..90b1c09 100644 --- a/src/includes/AVL.INC +++ b/src/includes/AVL.INC @@ -21,7 +21,7 @@ C NRMAX number of stored run cases C NTMAX number of stored time levels C C ICONX max number of control or design variable declaration lines per section (20) - PARAMETER (NVMAX=5550, + PARAMETER (NVMAX=5450, & NSMAX=400, & NFMAX=30, & NLMAX=500, @@ -54,7 +54,7 @@ C---- unit values, names, and namelengths COMMON /UN_C/UNCHL,UNCHM,UNCHT,UNCHF,UNCHS,UNCHV,UNCHA,UNCHI,UNCHD COMMON /UN_I/NUL ,NUM ,NUT ,NUF ,NUS ,NUV ,NUA ,NUI ,NUD C - CHARACTER*120 FRNDEF, FPRDEF, FEVDEF + CHARACTER*256 FILDEF, FRNDEF, FMSDEF, FPRDEF, FEVDEF CHARACTER*120 TITLE CHARACTER*40 STITLE, BTITLE, RTITLE CHARACTER*120 AFILES @@ -67,7 +67,9 @@ C CHARACTER*10 PARNAM CHARACTER*32 PARUNCH COMMON /CASE_C/ + & FILDEF, ! default configuration save file & FRNDEF, ! default run case save file + & FMSDEF, ! default mass distribution file & FPRDEF, ! default dimensional parameter file & FEVDEF, ! default eigenvalue save file & TITLE, ! configuration title @@ -91,6 +93,7 @@ C & LUOUT, ! logical unit for output dump file & LUSTD, ! logical unit for stability deriv. dump file & LUSYS, ! logical unit for dynamic system matrix dump file + & LUMAS, ! logical unit for mass file & IYSYM,IZSYM, ! y,z image symm. (0=no image, 1=image) & MATSYM, ! matrix symmetry flag & NVOR, ! number of horseshoe vortices @@ -110,8 +113,9 @@ C & IRUNT, ! target run case for time march initial state & ITRIM(NRMAX), ! trim type used for run case (if any) & NEIGEN(NRMAX), ! number of valid eigenmodes available for run case - & NEIGENDAT(NRMAX), ! number reference data eigenvalues - & CPSCL ! CP scaling factor for cp plot + & NEIGENDAT(NRMAX), ! number reference data eigenvalues + & CPSCL ! CP scaling factor for cp plot + C LOGICAL LGEO,LENC,LAIC,LSRD,LVEL,LSOL,LSEN, & LVISC,LMASS, @@ -211,13 +215,14 @@ C real(kind=avl_real) CHINGE_G real(kind=avl_real) DCL_A0, DCM_A0 real(kind=avl_real) DCL_U0, DCM_U0 - real(kind=avl_real) XNP + real(kind=avl_real) XNP + real(kind=avl_real) EXEC_TOL COMMON /CASE_R/ & VERSION, ! AVL version number & DTR, PI, ! 3.14159/180 , 3.14159 & YSYM, ZSYM, ! y- and z-locations of symmetry planes & ALFA, BETA, ! alpha, beta - & VINF(3), ! freestream velocies in body axes + & VINF(3), ! freestream velocity in body axes & VINF_A(3), ! d(Vinf)/d(alpha) & VINF_B(3), ! d(Vinf)/d(beta) & WROT(3), ! rotation rates in body axes @@ -241,7 +246,7 @@ C & CLFF_G(NGMAX),CYFF_G(NGMAX),CDFF_G(NGMAX), ! deriv wrt design & SPANEF, ! span efficiency & SPANEF_A, ! d(SPANEF)/d(alpha) - & SPANEF_U(NUMAX), ! d(SPANEF)/d(beta) + & SPANEF_U(NUMAX), ! d(SPANEF)/d(U,W)) & SPANEF_D(NDMAX), ! d(SPANEF)/d(control) & SPANEF_G(NGMAX), ! d(SPANEF)/d(design) & CDTOT, CLTOT, ! total CD,CL @@ -264,13 +269,14 @@ C & CDTOT_RX, CLTOT_RX, CYTOT_RX, CRTOT_RX, CMTOT_RX, CNTOT_RX, ! stab derives wrt beta & CDTOT_RY, CLTOT_RY, CYTOT_RY, CRTOT_RY, CMTOT_RY, CNTOT_RY, ! stab derives wrt beta & CDTOT_RZ, CLTOT_RZ, CYTOT_RZ, CRTOT_RZ, CMTOT_RZ, CNTOT_RZ, ! stab derives wrt beta - & CHINGE(NDMAX), - & CHINGE_U(NDMAX,NUMAX), - & CHINGE_D(NDMAX,NDMAX), - & CHINGE_G(NDMAX,NGMAX), + & CHINGE(NDMAX), !hinge moment for control + & CHINGE_U(NDMAX,NUMAX), ! sens wrt U,W + & CHINGE_D(NDMAX,NDMAX), ! sens wrt control + & CHINGE_G(NDMAX,NGMAX), ! sens wrt design & DCL_A0, DCM_A0, ! additional default CL_a, CM_a & DCL_U0, DCM_U0, ! additional default CL_u, CM_u - & XNP ! x neutral point + & XNP, ! x neutral point + & EXEC_TOL ! tolerance of exec trim solve COMPLEX EVAL, EVEC, EVALDAT COMMON /CASE_Z/ & EVAL(JEMAX,NRMAX), ! mode eigenvalue @@ -311,11 +317,12 @@ C & AINER(3,3) ! apparent inertia/rho | from geometry - LOGICAL LFWAKE, LFALBE, LFLOAD + LOGICAL LFWAKE, LFALBE, LFLOAD, LRANGE COMMON /SURF_L/ & LFWAKE(NFMAX), ! T if surface sheds a wake & LFALBE(NFMAX), ! T if surface is to see freestream alpha,beta - & LFLOAD(NFMAX) ! T if surface contributes to overall loads + & LFLOAD(NFMAX), ! T if surface contributes to overall loads + & LRANGE(NFMAX) ! T if surface determined using full airfoil range COMMON /SURF_I/ & NJ(NFMAX), ! number of elements along span in surface @@ -323,7 +330,13 @@ C & IFRST(NFMAX), ! index of first element in surface & JFRST(NFMAX), ! index of first strip in surface & IMAGS(NFMAX), ! indicates whether surface is a YDUPlicated one - & LSCOMP(NFMAX) ! logical surface component index +cc#ifdef USE_CPOML + & LSCOMP(NFMAX), ! logical surface component index + & ICNTFRST(NFMAX), ! index of first section counter in surface + & NCNTSEC(NFMAX) ! number of section counters in surface +cc#else +cc & LSCOMP(NFMAX) ! logical surface component index +cc#endif REAL(kind=avl_real) CDSURF,CLSURF REAL(kind=avl_real) CXSURF,CYSURF,CZSURF @@ -344,19 +357,20 @@ C REAL(kind=avl_real) SSURF, CAVESURF REAL(kind=avl_real) VBODY COMMON /SURF_R/ +C Surface forces and moments referenced to SREF,CREF,BREF about XREF,YREF,ZREF & CDSURF(NFMAX),CLSURF(NFMAX), ! surface CD,CL & CXSURF(NFMAX),CYSURF(NFMAX),CZSURF(NFMAX), ! surface Cx,Cy,Cz & CRSURF(NFMAX),CNSURF(NFMAX),CMSURF(NFMAX), ! surface Cl,Cm,Cn & CDVSURF(NFMAX), ! surface viscous CD - & CDS_A(NFMAX), CLS_A(NFMAX), + & CDS_A(NFMAX), CLS_A(NFMAX), ! alpha sens. & CDS_U(NFMAX,NUMAX), CLS_U(NFMAX,NUMAX), - & CXS_U(NFMAX,NUMAX), CYS_U(NFMAX,NUMAX), CZS_U(NFMAX,NUMAX), + & CXS_U(NFMAX,NUMAX), CYS_U(NFMAX,NUMAX), CZS_U(NFMAX,NUMAX), ! velocity and rotation sens. & CRS_U(NFMAX,NUMAX), CNS_U(NFMAX,NUMAX), CMS_U(NFMAX,NUMAX), & CDS_D(NFMAX,NDMAX), CLS_D(NFMAX,NDMAX), - & CXS_D(NFMAX,NDMAX), CYS_D(NFMAX,NDMAX), CZS_D(NFMAX,NDMAX), + & CXS_D(NFMAX,NDMAX), CYS_D(NFMAX,NDMAX), CZS_D(NFMAX,NDMAX), ! control sens. & CRS_D(NFMAX,NDMAX), CNS_D(NFMAX,NDMAX), CMS_D(NFMAX,NDMAX), & CDS_G(NFMAX,NGMAX), CLS_G(NFMAX,NGMAX), - & CXS_G(NFMAX,NGMAX), CYS_G(NFMAX,NGMAX), CZS_G(NFMAX,NGMAX), + & CXS_G(NFMAX,NGMAX), CYS_G(NFMAX,NGMAX), CZS_G(NFMAX,NGMAX), ! design (incidence) sens. & CRS_G(NFMAX,NGMAX), CNS_G(NFMAX,NGMAX), CMS_G(NFMAX,NGMAX), & CF_SRF(3,NFMAX), CM_SRF(3,NFMAX), & CL_SRF(NFMAX), CD_SRF(NFMAX), CMLE_SRF(NFMAX), @@ -410,12 +424,17 @@ c REAL(kind=avl_real) SSPACE REAL(kind=avl_real) SSPACES REAL(kind=avl_real) XYZLES + REAL(kind=avl_real) XLASEC + REAL(kind=avl_real) ZLASEC + REAL(kind=avl_real) XUASEC + REAL(kind=avl_real) ZUASEC REAL(kind=avl_real) CHORDS REAL(kind=avl_real) AINCS REAL(kind=avl_real) XASEC REAL(kind=avl_real) SASEC REAL(kind=avl_real) TASEC REAL(kind=avl_real) CLCDSEC + REAL(kind=avl_real) CLCDSRF REAL(kind=avl_real) CLAF REAL(kind=avl_real) XHINGED REAL(kind=avl_real) VHINGED @@ -431,12 +450,17 @@ c & SSPACE(NFMAX), ! spanwise spacing & SSPACES(NSMAX, NFMAX), ! spanwise spacing vector & XYZLES(3, NSMAX, NFMAX), ! leading edge cordinate vector + & XLASEC(IBX, NSMAX, NFMAX), + & ZLASEC(IBX, NSMAX, NFMAX), + & XUASEC(IBX, NSMAX, NFMAX), + & ZUASEC(IBX, NSMAX, NFMAX), & CHORDS(NSMAX, NFMAX), ! chord length vectorm & AINCS(NSMAX, NFMAX), ! incidence angle vector & XASEC(IBX, NSMAX, NFMAX), ! the x coordinate aifoil section & SASEC(IBX, NSMAX, NFMAX), ! slope of camber line at xasec & TASEC(IBX, NSMAX, NFMAX), ! thickness at xasec & CLCDSEC(6, NSMAX, NFMAX), + & CLCDSRF(6, NFMAX), & CLAF(NSMAX, NFMAX), & XHINGED(ICONX, NSMAX, NFMAX), ! hinge location & VHINGED(3, ICONX, NSMAX, NFMAX), ! hinge vector @@ -450,7 +474,8 @@ C !!--- end added variables for python geometry minipulation --- COMMON /STRP_I/ & NSURFS(NSMAX), ! index of surface which contains this strip & IJFRST(NSMAX), ! index of first element in strip - & NVSTRP(NSMAX) ! number of elements in strip + & NVSTRP(NSMAX), ! number of elements in strip + & ICNTSEC(NSMAX) ! section counters for surface LOGICAL LSTRIPOFF,LVISCSTRP,LJ1SECT,LJ2SECT COMMON /STRP_L/ @@ -503,31 +528,43 @@ C !!--- end added variables for python geometry minipulation --- & SAXFR, ! x/c of spanwise axis for Vperp def & ESS(3,NSMAX), ! spanwise unit vector for Vperp def & ENSY(NSMAX), ENSZ(NSMAX), ! strip normal vector in Trefftz-Pln - & XSREF(NSMAX),YSREF(NSMAX),ZSREF(NSMAX), + & XSREF(NSMAX),YSREF(NSMAX),ZSREF(NSMAX), ! strip reference point & AINC(NSMAX), ! strip's incidence twist angle & AINC_G(NSMAX,NGMAX), ! dAINC/dG - & CDSTRP(NSMAX), CLSTRP(NSMAX), ! strip forces - & CXSTRP(NSMAX), CYSTRP(NSMAX), CZSTRP(NSMAX), ! strip forces - & CRSTRP(NSMAX), CNSTRP(NSMAX), CMSTRP(NSMAX), ! strip moments +C Strip forces and moments + & CDSTRP(NSMAX), CLSTRP(NSMAX), ! strip forces in stability axes + & CXSTRP(NSMAX), CYSTRP(NSMAX), CZSTRP(NSMAX), ! strip forces in body axes + & CRSTRP(NSMAX), CNSTRP(NSMAX), CMSTRP(NSMAX), ! strip moments in body axes about XREF,YREF,ZREF +C & CDST_A(NSMAX), CLST_A(NSMAX), ! alpha sens. & CDST_U(NSMAX,NUMAX), CLST_U(NSMAX,NUMAX), - & CXST_U(NSMAX,NUMAX), CYST_U(NSMAX,NUMAX), CZST_U(NSMAX,NUMAX), + & CXST_U(NSMAX,NUMAX), CYST_U(NSMAX,NUMAX), CZST_U(NSMAX,NUMAX), ! freestream velocity and rotation sens. & CRST_U(NSMAX,NUMAX), CNST_U(NSMAX,NUMAX), CMST_U(NSMAX,NUMAX), & CDST_D(NSMAX,NDMAX), CLST_D(NSMAX,NDMAX), - & CXST_D(NSMAX,NDMAX), CYST_D(NSMAX,NDMAX), CZST_D(NSMAX,NDMAX), + & CXST_D(NSMAX,NDMAX), CYST_D(NSMAX,NDMAX), CZST_D(NSMAX,NDMAX), ! control sens. & CRST_D(NSMAX,NDMAX), CNST_D(NSMAX,NDMAX), CMST_D(NSMAX,NDMAX), & CDST_G(NSMAX,NGMAX), CLST_G(NSMAX,NGMAX), - & CXST_G(NSMAX,NGMAX), CYST_G(NSMAX,NGMAX), CZST_G(NSMAX,NGMAX), + & CXST_G(NSMAX,NGMAX), CYST_G(NSMAX,NGMAX), CZST_G(NSMAX,NGMAX), ! design (incidence) sens. & CRST_G(NSMAX,NGMAX), CNST_G(NSMAX,NGMAX), CMST_G(NSMAX,NGMAX), - & CF_STRP(3,NSMAX), CM_STRP(3,NSMAX), ! strip forces, norm local - & CNRMSTRP(NSMAX), CAXLSTRP(NSMAX), - & CD_LSTRP(NSMAX), CL_LSTRP(NSMAX), - & CDV_LSTRP(NSMAX), CLTSTRP(NSMAX), CLASTRP(NSMAX), - & CMC4(NSMAX), CMLE(NSMAX), - & CNC(NSMAX), DWWAKE(NSMAX), - & CNC_U(NSMAX,NUMAX), +C + & CF_STRP(3,NSMAX), CM_STRP(3,NSMAX), ! strip forces in body axes referenced to strip area and 1/4 chord + & CNRMSTRP(NSMAX), CAXLSTRP(NSMAX), ! strip forces in local dihedral plane (normal and axial forces) + & CD_LSTRP(NSMAX), CL_LSTRP(NSMAX), ! strip forces in local dihedral plane (lift and drag forces) + & CDV_LSTRP(NSMAX), ! strip viscous drag in stability axes + & CLTSTRP(NSMAX), CLASTRP(NSMAX), ! strip CL referenced to Vperp, CL referenced to total local velocity + & CMC4(NSMAX), CMLE(NSMAX), ! strip pitching moment about c/4 and pitching moment about LE vector + & CNC(NSMAX), DWWAKE(NSMAX), ! strip spanloading and downwash in Trefftz plane + & CNC_U(NSMAX,NUMAX), ! spanloading sens. & CNC_D(NSMAX,NDMAX), & CNC_G(NSMAX,NGMAX) +C +C +cc#ifdef USE_CPOML + REAL(kind=avl_real) AINC1 + REAL(kind=avl_real) AINC2 + COMMON /STRP_S/ + & AINC1(NSMAX), AINC2(NSMAX) ! left/right strip incidence twist angle +cc#endif C LOGICAL LVNC, LVALBE COMMON /VRTX_L/ @@ -577,6 +614,7 @@ C REAL(kind=avl_real) RHS REAL(kind=avl_real) RES REAL(kind=avl_real) RES_D + REAL(kind=avl_real) CPT COMMON /VRTX_R/ & RV1(3,NVMAX), ! h.v. vortex left points & RV2(3,NVMAX), ! h.v. vortex right points @@ -588,36 +626,52 @@ C & CHORDV(NVMAX), ! chord of element-containing strip & SLOPEV(NVMAX), ! camber slopes at h.v. bound leg & SLOPEC(NVMAX), ! camber slopes at c.p. - & DCONTROL(NVMAX,NDMAX), ! d(normal angle)/dCONTROL - & VHINGE(3,NSMAX,NDMAX), ! hinge vector for CONTROL rot. of normal - & PHINGE(3,NSMAX,NDMAX), ! point on hingeline for hinge moment calculation - & VREFL(NSMAX,NDMAX), !sign applied to hinge vec. of refl. surface - & ENC(3,NVMAX), ! c.p. normal vector - & ENV(3,NVMAX), ! h.v. normal vector - & ENC_D(3,NVMAX,NDMAX), ! sensitivities + & DCONTROL(NVMAX,NDMAX), ! d(normal angle)/dCONTROL + & VHINGE(3,NSMAX,NDMAX), ! hinge vector for CONTROL rot. of normal + & PHINGE(3,NSMAX,NDMAX), ! point on hingeline for hinge moment calculation + & VREFL(NSMAX,NDMAX), ! sign applied to hinge vec. of refl. surface + & ENC(3,NVMAX), ! control point normal vector + & ENV(3,NVMAX), ! horseshoe vortex normal vector + & ENC_D(3,NVMAX,NDMAX), ! control point normal vector sensitivities & ENC_G(3,NVMAX,NGMAX), - & ENV_D(3,NVMAX,NDMAX), + & ENV_D(3,NVMAX,NDMAX), ! horseshoe vortex normal vector sensitivities & ENV_G(3,NVMAX,NGMAX), - & DCP(NVMAX), ! delta(Cp) on vortex element + & DCP(NVMAX), ! delta(Cp) on vortex element & DCP_U(NVMAX,NUMAX), & DCP_D(NVMAX,NDMAX), & DCP_G(NVMAX,NGMAX), - & GAM(NVMAX), ! circulation of h.v. vortex + & GAM(NVMAX), ! circulation of horseshoe vortex & GAM_U_0(NVMAX,NUMAX), & GAM_U_D(NVMAX,NUMAX,NDMAX), & GAM_U_G(NVMAX,NUMAX,NGMAX), - & GAM_U(NVMAX,NUMAX), + & GAM_U(NVMAX,NUMAX), ! circulation sensitivities & GAM_D(NVMAX,NDMAX), & GAM_G(NVMAX,NGMAX), - & SRC(NLMAX), ! source strength of source+doublet line elem - & DBL(3,NLMAX), ! doublet strength of source+doublet line elem - & SRC_U(NLMAX,NUMAX), + & SRC(NLMAX), ! source strength of source+doublet line elem + & DBL(3,NLMAX), ! doublet strength of source+doublet line elem + & SRC_U(NLMAX,NUMAX), ! sensitivities & DBL_U(3,NLMAX,NUMAX), & WCSRD_U(3,NVMAX,NUMAX), +cc#ifdef USE_CPOML & WVSRD_U(3,NVMAX,NUMAX), - & RHS(NVMAX), ! RHS (flow normal to surface) for linear system - & RES(NVMAX), ! residual for linear system res = Ax-b - & RES_D(NVMAX,NDMAX) ! residual for gamma_d lin. sys + & CPT(NVMAX), ! thickness-based Cp at panel midpoint +cc#else + & RHS(NVMAX), ! RHS (flow normal to surface) for linear system + & RES(NVMAX), ! residual for linear system res = Ax-b + & RES_D(NVMAX,NDMAX) ! residual for gamma_d lin. sys + +cc#ifdef USE_CPOML + REAL(kind=avl_real) XYN1 + REAL(kind=avl_real) XYN2 + REAL(kind=avl_real) ZLON1 + REAL(kind=avl_real) ZLON2 + REAL(kind=avl_real) ZUPN1 + REAL(kind=avl_real) ZUPN2 + COMMON /VRTX_S/ + & XYN1(2,NVMAX), XYN2(2,NVMAX), ! left/right aft-node of element + & ZLON1(NVMAX), ZLON2(NVMAX), ! left/right lower z-coord of aft-node + & ZUPN1(NVMAX), ZUPN2(NVMAX) ! left/right upper z-coord of aft-node +cc#endif C COMMON /BODY_I/ & NL(NBMAX), ! number of source-line nodes in body @@ -631,9 +685,9 @@ C REAL(kind=avl_real) CXBDY,CYBDY,CZBDY REAL(kind=avl_real) CRBDY,CNBDY,CMBDY COMMON /BODY_R/ - & ELBDY(NBMAX), ! body length - & SRFBDY(NBMAX), ! body surface area - & VOLBDY(NBMAX), ! body volume + & ELBDY(NBMAX), ! body length + & SRFBDY(NBMAX), ! body surface area + & VOLBDY(NBMAX), ! body volume & DCPB(3,NLMAX), ! dCp loading on sources & CDBDY(NBMAX),CLBDY(NBMAX), ! body CD,CL & CXBDY(NBMAX),CYBDY(NBMAX),CZBDY(NBMAX), ! body CX,CY,CZ @@ -689,7 +743,8 @@ C & LPLTSURF(NFMAX), & LPLTBODY(NBMAX), & LSVMOV - COMMON /PLOT_I/ NTRI, IMARKSURF, NAXANN(3), IRCOLOR(NRMAX) + COMMON /PLOT_I/ NTRI, IMARKSURF, ICOLORSURF, NAXANN(3), + & IRCOLOR(NRMAX) REAL(kind=avl_real) TRI REAL(kind=avl_real) AXMIN,AXMAX,AXDEL,AXSPAN diff --git a/src/includes/AVLPLT.INC b/src/includes/AVLPLT.INC index a88c2d5..f35586f 100644 --- a/src/includes/AVLPLT.INC +++ b/src/includes/AVLPLT.INC @@ -12,6 +12,7 @@ C & XMARG,YMARG, & XPAGE,YPAGE, & XWIND,YWIND, + & UNT(3,3), ORG(3,3), & PMARG,CPFAC, HNFAC, ENFAC, & AZIMOB,ELEVOB,TILTOB,ROBINV C diff --git a/src/unitset.f b/src/unitset.f new file mode 100644 index 0000000..5a6794c --- /dev/null +++ b/src/unitset.f @@ -0,0 +1,486 @@ +C*********************************************************************** +C Module: amass.f +C +C Copyright (C) 2002 Mark Drela, Harold Youngren +C +C This program is free software; you can redistribute it and/or modify +C it under the terms of the GNU General Public License as published by +C the Free Software Foundation; either version 2 of the License, or +C (at your option) any later version. +C +C This program is distributed in the hope that it will be useful, +C but WITHOUT ANY WARRANTY; without even the implied warranty of +C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +C GNU General Public License for more details. +C +C You should have received a copy of the GNU General Public License +C along with this program; if not, write to the Free Software +C Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. +C*********************************************************************** + + SUBROUTINE MASINI + INCLUDE 'AVL.INC' +C--------------------------------------------------- +C Initializes default mass, inertia tensor +C--------------------------------------------------- + RMASS0 = 1.0 +C + DO K = 1, 3 + RINER0(K,1) = 0. + RINER0(K,2) = 0. + RINER0(K,3) = 0. + RINER0(K,K) = 1.0 +C + AMASS(K,1) = 0. + AMASS(K,2) = 0. + AMASS(K,3) = 0. +C + AINER(K,1) = 0. + AINER(K,2) = 0. + AINER(K,3) = 0. + ENDDO +C + XYZMASS0(1) = 0. + XYZMASS0(2) = 0. + XYZMASS0(3) = 0. +C + LMASS = .FALSE. +C + RETURN + END ! MASINI + + + + SUBROUTINE MASGET(LU,FNAME,ERROR) +C-------------------------------------------------- +C Reads mass distributions file, +C computes default mass, inertia tensor +C-------------------------------------------------- + INCLUDE 'AVL.INC' + CHARACTER*(*) FNAME +C + REAL mass, mi + REAL Ixx , Iyy , Izz , Ixy , Ixz , Iyz + REAL Ixxi, Iyyi, Izzi, Ixyi, Ixzi, Iyzi +C + CHARACTER*80 CASE + CHARACTER*256 LINE, LINEU +C + CHARACTER*32 UNCHGEE, UNCHRHO, UNCH +C + LOGICAL ERROR + REAL RINP(10), FAC(10), ADD(10) +C + 1000 FORMAT(A) +C + OPEN(LU,FILE=FNAME,STATUS='OLD',ERR=90) +C + sum_m = 0. + sum_mx = 0. + sum_my = 0. + sum_mz = 0. + sum_mxx = 0. + sum_myy = 0. + sum_mzz = 0. + sum_mxy = 0. + sum_mxz = 0. + sum_myz = 0. + sum_mzz = 0. + sum_ixx = 0. + sum_iyy = 0. + sum_izz = 0. + sum_ixy = 0. + sum_iyz = 0. + sum_ixz = 0. +C + xcg = 0. + ycg = 0. + zcg = 0. +C +C---- default multipliers and adders + DO K = 1, 10 + FAC(K) = 1.0 + ADD(K) = 0.0 + ENDDO +C + UNITL = 1. + UNITM = 1. + UNITT = 1. + UNCHL = 'Lunit' + UNCHM = 'Munit' + UNCHT = 'Tunit' + NUL = 5 + NUM = 5 + NUT = 5 +C + GEE0 = 1. + RHO0 = 1. + UNCHGEE = 'Lunit/Tunit^2' + UNCHRHO = 'Munit/Lunit^3' + NUGEE = 13 + NURHO = 13 +C + ILINE = 0 + +C---- search all lines of file for data values +C============================================= + 10 CONTINUE + READ(LU,1000,END=50) LINE + ILINE = ILINE + 1 +C + IF(INDEX('#!',LINE(1:1)).NE.0) GO TO 10 +C + IF(INDEX('*',LINE(1:1)).NE.0) THEN +C------- read multiplier line + KINP = 10 + CALL GETFLT(LINE(2:80),RINP,KINP,ERROR) + IF(ERROR) GO TO 40 +C + DO K = 1, KINP + FAC(K) = RINP(K) + ENDDO + GO TO 10 + ENDIF +C + IF(INDEX('+',LINE(1:1)).NE.0) THEN +C------- read adder line + KINP = 10 + CALL GETFLT(LINE(2:80),RINP,KINP,ERROR) + IF(ERROR) GO TO 40 +C + DO K = 1, KINP + ADD(K) = RINP(K) + ENDDO + GO TO 10 + ENDIF +C +C +C------ don't look for special parameters if there's no "=" character + KEQ = INDEX(LINE,'=') + IF(KEQ.LE.1) GO TO 20 +C +C------------------------------------------ +C------ look for parameter keywords in front of "=" character +C + K = INDEX(LINE(1:KEQ-1),'Lunit') + IF(K .NE. 0) THEN + LINEU = LINE(KEQ+1:256) + CALL STRIP(LINEU,NLN) + KB = INDEX(LINEU,' ') + READ(LINEU(1:KB),*,ERR=40) UNITL + UNCHL = LINEU(KB+1:256) + CALL STRIP(UNCHL,NUL) + NUL = MAX(NUL,1) + GO TO 10 + ENDIF +C + K = INDEX(LINE(1:KEQ-1),'Munit') + IF(K .NE. 0) THEN + LINEU = LINE(KEQ+1:256) + CALL STRIP(LINEU,NLN) + KB = INDEX(LINEU,' ') + READ(LINEU(1:KB),*,ERR=40) UNITM + UNCHM = LINEU(KB+1:256) + CALL STRIP(UNCHM,NUM) + NUM = MAX(NUM,1) + GO TO 10 + ENDIF +C + K = INDEX(LINE(1:KEQ-1),'Tunit') + IF(K .NE. 0) THEN + LINEU = LINE(KEQ+1:256) + CALL STRIP(LINEU,NLN) + KB = INDEX(LINEU,' ') + READ(LINEU(1:KB),*,ERR=40) UNITT + UNCHT = LINEU(KB+1:256) + CALL STRIP(UNCHT,NUT) + NUT = MAX(NUT,1) + GO TO 10 + ENDIF +C +C + K = INDEX(LINE(1:KEQ-1),'g') + IF(K .NE. 0) THEN + LINEU = LINE(KEQ+1:256) + CALL STRIP(LINEU,NLN) + KB = INDEX(LINEU,' ') + READ(LINEU(1:KB),*,ERR=40) GEE0 + UNCHGEE = LINEU(KB+1:256) + CALL STRIP(UNCHGEE,NUGEE) + NUGEE = MAX(NUGEE,1) + GO TO 10 + ENDIF +C + K = INDEX(LINE(1:KEQ-1),'rho') + IF(K .NE. 0) THEN + LINEU = LINE(KEQ+1:256) + CALL STRIP(LINEU,NLN) + KB = INDEX(LINEU,' ') + READ(LINEU(1:KB),*,ERR=40) RHO0 + UNCHRHO = LINEU(KB+1:256) + CALL STRIP(UNCHRHO,NURHO) + NURHO = MAX(NURHO,1) + GO TO 10 + ENDIF +C +C------------------------------ + 20 CONTINUE +C + DO I=1, 10 + RINP(I) = 0. + ENDDO +C + KINP = 10 + CALL GETFLT(LINE,RINP,KINP,ERROR) + IF(ERROR) GO TO 40 +C + mi = FAC(1)*RINP(1) + ADD(1) + xi = FAC(2)*RINP(2) + ADD(2) + yi = FAC(3)*RINP(3) + ADD(3) + zi = FAC(4)*RINP(4) + ADD(4) + Ixxi = FAC(5)*RINP(5) + ADD(5) + Iyyi = FAC(6)*RINP(6) + ADD(6) + Izzi = FAC(7)*RINP(7) + ADD(7) + Ixyi = FAC(8)*RINP(8) + ADD(8) + Ixzi = FAC(9)*RINP(9) + ADD(9) + Iyzi = FAC(10)*RINP(10) + ADD(10) +C + sum_m = sum_m + mi + sum_mx = sum_mx + mi*xi + sum_my = sum_my + mi*yi + sum_mz = sum_mz + mi*zi + sum_mxx = sum_mxx + mi*xi*xi + sum_myy = sum_myy + mi*yi*yi + sum_mzz = sum_mzz + mi*zi*zi + sum_mxy = sum_mxy + mi*xi*yi + sum_mxz = sum_mxz + mi*xi*zi + sum_myz = sum_myz + mi*yi*zi +c + sum_ixx = sum_ixx + Ixxi + sum_iyy = sum_iyy + Iyyi + sum_izz = sum_izz + Izzi + sum_ixy = sum_ixy + Ixyi + sum_ixz = sum_ixz + Ixzi + sum_iyz = sum_iyz + Iyzi +C + GO TO 10 +C +C----------------------------------- + 40 WRITE(*,1500) ILINE, LINE(1:80) + 1500 FORMAT(' * Bad data line',I5,' ignored: ',A) + GO TO 10 +C +C============================================= +C + 50 CONTINUE + CLOSE(LU) +C + mass = sum_m +C + xcg = sum_mx/sum_m + ycg = sum_my/sum_m + zcg = sum_mz/sum_m +C + Ixx = sum_ixx + (sum_myy + sum_mzz) - sum_m*(ycg**2 + zcg**2) + Iyy = sum_iyy + (sum_mzz + sum_mxx) - sum_m*(zcg**2 + xcg**2) + Izz = sum_izz + (sum_mxx + sum_myy) - sum_m*(xcg**2 + ycg**2) + Ixy = sum_ixy + sum_mxy - sum_m* xcg*ycg + Ixz = sum_ixz + sum_mxz - sum_m* xcg*zcg + Iyz = sum_iyz + sum_myz - sum_m* ycg*zcg +C + RMASS0 = mass * UNITM +C + RINER0(1,1) = Ixx * UNITM*UNITL**2 + RINER0(1,2) = -Ixy * UNITM*UNITL**2 + RINER0(1,3) = -Ixz * UNITM*UNITL**2 +C + RINER0(2,1) = -Ixy * UNITM*UNITL**2 + RINER0(2,2) = Iyy * UNITM*UNITL**2 + RINER0(2,3) = -Iyz * UNITM*UNITL**2 +C + RINER0(3,1) = -Ixz * UNITM*UNITL**2 + RINER0(3,2) = -Iyz * UNITM*UNITL**2 + RINER0(3,3) = Izz * UNITM*UNITL**2 +C + XYZMASS0(1) = xcg * UNITL + XYZMASS0(2) = ycg * UNITL + XYZMASS0(3) = zcg * UNITL +C +C---------------------------------------- +C---- set derived unit values and names + CALL UNITSET +C +C---- set parameter unit names + CALL PARNSET +C + ERROR = .FALSE. + LMASS = .TRUE. + RETURN +C + 90 CONTINUE + CALL STRIP(FNAME,NFN) + WRITE(*,9000) FNAME(1:NFN) + 9000 FORMAT(/' Mass file ',A,' open error') + ERROR = .TRUE. + LMASS = .FALSE. + RETURN + END ! MASGET + + + + SUBROUTINE MASPUT(IR1,IR2) + INCLUDE 'AVL.INC' +C------------------------------------------- +C Stores default mass, inertias +C in run case parameter array +C------------------------------------------- +C + DO IR = IR1, IR2 + PARVAL(IPMASS,IR) = RMASS0 + PARVAL(IPIXX,IR) = RINER0(1,1) + PARVAL(IPIYY,IR) = RINER0(2,2) + PARVAL(IPIZZ,IR) = RINER0(3,3) + PARVAL(IPIXY,IR) = RINER0(1,2) + PARVAL(IPIYZ,IR) = RINER0(2,3) + PARVAL(IPIZX,IR) = RINER0(3,1) + + PARVAL(IPGEE,IR) = GEE0 + PARVAL(IPRHO,IR) = RHO0 +C + PARVAL(IPXCG,IR) = XYZMASS0(1)/UNITL + PARVAL(IPYCG,IR) = XYZMASS0(2)/UNITL + PARVAL(IPZCG,IR) = XYZMASS0(3)/UNITL + ENDDO +C + RETURN + END ! MASPUT + + + + SUBROUTINE MASSHO(LU) + INCLUDE 'AVL.INC' +C + WRITE(LU,*) + IF(UNCHM(1:NUM).NE.'Munit') + &WRITE(LU,1250) 'Mass = ', RMASS0/UNITM, 'Munit' + WRITE(LU,1250) 'Mass = ', RMASS0, UNCHM(1:NUM) +C + WRITE(LU,*) + WRITE(LU,1260) 'Ref. x,y,z = ',(XYZREF0(K) , K=1,3),'Lunit' + IF(UNCHL(1:NUL).NE.'Lunit') + &WRITE(LU,1260) 'C.G. x,y,z = ',(XYZMASS0(K)/UNITL, K=1,3),'Lunit' + WRITE(LU,1260) 'C.G. x,y,z = ',(XYZMASS0(K), K=1,3), UNCHL(1:NUL) +C + WRITE(LU,*) + WRITE(LU,1271) 'Ixx -Ixy -Ixz | ', + & RINER0(1,1),RINER0(1,2),RINER0(1,3),'|' + WRITE(LU,1272) ' Iyy -Iyz = | ', + & RINER0(2,2),RINER0(2,3),'|', UNCHI(1:NUI) + WRITE(LU,1273) ' Izz | ', + & RINER0(3,3),'|' +C + 1250 FORMAT(1X, A, G12.4,2X,A) + 1260 FORMAT(1X, A, 3G12.4,2X,A) + 1271 FORMAT(1X, A, 3G12.4, 2X, A, 2X, A) + 1272 FORMAT(1X, A, 12X, 2G12.4, 2X, A, 2X, A) + 1273 FORMAT(1X, A, 24X, G12.4, 2X, A, 2X, A) +C + RETURN + END ! MASSHO + + + SUBROUTINE APPSHO(LU,RHO) + INCLUDE 'AVL.INC' +C + WRITE(LU,*) 'Apparent mass, inertia' + WRITE(LU,*) + WRITE(LU,1271) 'mxx mxy mxz | ', + & AMASS(1,1)*RHO,AMASS(1,2)*RHO,AMASS(1,3)*RHO,'|' + WRITE(LU,1272) ' myy myz = | ', + & AMASS(2,2)*RHO,AMASS(2,3)*RHO,'|', UNCHM(1:NUM) + WRITE(LU,1273) ' mzz | ', + & AMASS(3,3)*RHO,'|' + WRITE(LU,*) + WRITE(LU,1271) 'Ixx -Ixy -Ixz | ', + & AINER(1,1)*RHO,AINER(1,2)*RHO,AINER(1,3)*RHO,'|' + WRITE(LU,1272) ' Iyy -Iyz = | ', + & AINER(2,2)*RHO,AINER(2,3)*RHO,'|', UNCHI(1:NUI) + WRITE(LU,1273) ' Izz | ', + & AINER(3,3)*RHO,'|' +C + 1271 FORMAT(1X, A, 3G12.4, 2X, A, 2X, A) + 1272 FORMAT(1X, A, 12X, 2G12.4, 2X, A, 2X, A) + 1273 FORMAT(1X, A, 24X, G12.4, 2X, A, 2X, A) +C + RETURN + END ! APPSHO + + + + SUBROUTINE APPGET +C---------------------------------------------------------------- +C Calculates apparent mass and inertia for unit air density +C---------------------------------------------------------------- + INCLUDE 'AVL.INC' + REAL UC(3), US(3), UN(3), RM(3), RXUN(3) +C + DO K = 1, 3 + DO L = 1, 3 + AMASS(K,L) = 0. + AINER(K,L) = 0. + ENDDO + ENDDO +C + DO 100 J = 1, NSTRIP + CR = CHORD(J) + SR = CHORD(J)*WSTRIP(J) +C +C------ strip normal unit vector UN + UN(1) = 0.0 + UN(2) = ENSY(J) + UN(3) = ENSZ(J) +C +C------ spanwise unit vector US (along midchord) + US(1) = RLE2(1,J) - RLE1(1,J) + 0.5*(CHORD2(J) - CHORD1(J)) + US(2) = RLE2(2,J) - RLE1(2,J) + US(3) = RLE2(3,J) - RLE1(3,J) + UMAG = SQRT(US(1)*US(1) + US(2)*US(2) + US(3)*US(3)) + IF(UMAG.GT.0.0) THEN + US(1) = US(1)/UMAG + US(2) = US(2)/UMAG + US(3) = US(3)/UMAG + ENDIF +C +C------ perpendicular-chord unit vector UC + CALL CROSS(US,UN,UC) +C +C------ use the strip 1/2 chord location for moment arm + RM(1) = RLE(1,J) + 0.5*CR + RM(2) = RLE(2,J) + RM(3) = RLE(3,J) +C + CALL CROSS(RM,UN,RXUN) +C +C------ perpendicular chord + CPERP = CR * (US(2)*UN(3) - US(3)*UN(2)) +C +C------ apparent mass and spanwise-axis apparent inertia of strip + APPM = SR * 0.25*PI*CPERP + APPI = SR * 0.25*PI*CPERP**3 / 64.0 +C +C------ add apparent mass contribution + DO K = 1, 3 + DO L = 1, 3 + AMASS(K,L) = AMASS(K,L) + APPM* UN(K)* UN(L) * UNITL**3 + AINER(K,L) = AINER(K,L) + APPM*RXUN(K)*RXUN(L) * UNITL**5 + & + APPI* US(K)* US(L) * UNITL**5 + ENDDO + ENDDO +C + 100 CONTINUE +C + RETURN + END ! APPGET + + + diff --git a/tests/test_analysis.py b/tests/test_analysis.py index 7b91998..0e99b64 100644 --- a/tests/test_analysis.py +++ b/tests/test_analysis.py @@ -101,7 +101,8 @@ def test_constrained_cl_sweep(self): cl_arr = np.arange(0.6, 1.7, 0.1) for idx_cl, cl in enumerate(cl_arr): self.avl_solver.add_trim_condition("CL", cl) - self.avl_solver.execute_run() + # the tight tolerance here helps catch small issues with the newton solver + self.avl_solver.execute_run(tol=1e-12) run_data = self.avl_solver.get_case_total_data() np.testing.assert_allclose( @@ -130,9 +131,10 @@ def test_coefs(self): self.avl_solver.execute_run() coef_data = self.avl_solver.get_case_total_data() print(coef_data) + # the values are wonky here because of an unrealistic CDCL curve np.testing.assert_allclose(coef_data["CL"], 0.636031170179549, rtol=1e-8) - np.testing.assert_allclose(coef_data["CD"], 0.022842500250874982, rtol=1e-8) - np.testing.assert_allclose(coef_data["CM"], 0.00139609360952168, rtol=1e-8) + np.testing.assert_allclose(coef_data["CD"], 3.6953247032454204, rtol=1e-8) + np.testing.assert_allclose(coef_data["CM"], -0.5736410313952236, rtol=1e-8) class TestHingeMom(unittest.TestCase): def setUp(self): diff --git a/tests/test_io.py b/tests/test_io.py index c133dc1..060dc6c 100644 --- a/tests/test_io.py +++ b/tests/test_io.py @@ -46,7 +46,7 @@ def test_write_geom(self): avl_solver.write_geom_file(geom_output_file) baseline_data = avl_solver.get_surface_params() baseline_data_body = avl_solver.get_body_params() - # import pdb; pdb.set_trace() + del avl_solver avl_solver = AVLSolver(geo_file=geom_output_file) new_data = avl_solver.get_surface_params() @@ -114,8 +114,7 @@ def setUp(self): self.avl_solver = AVLSolver(geo_file=geom_mod_file, mass_file=mass_file) def test_get_scalar(self): - avl_version = 3.35 - + avl_version = 3.40 version = self.avl_solver.get_avl_fort_arr("CASE_R", "VERSION") self.assertEqual(version, avl_version)