.NET / Rotor源码分析5 - 开始使用WinDbg+SOS调试,sscoree.dll,加载SOS并设置JIT断点

本文详细介绍了如何使用WinDbg+SOS调试.NET/Rotor托管代码的过程,从准备调试环境到跟踪CLIX启动,再到加载SOS扩展,最终设置托管代码断点。

准备工作

在经过一番准备之后,现在我们可以开始正式使用WinDbg+SOS来调试托管代码了。如果你没有看过前两篇文章,那么请先阅读这两篇文章以对WinDbg+SOS有一个大致的了解。这两篇文章的链接在这里:

.NET Rotor源码研究4 – 修改Rotor使其发送CLR Notificationhttp://blog.youkuaiyun.com/ATField/archive/2007/05/21/1618535.aspx

.NET Rotor源码研究3 - 调试Rotor托管代码的利器:WinDbgSOShttp://blog.youkuaiyun.com/ATField/archive/2007/05/12/1606151.aspx

除此之外,还需要准备一个小程序来进行调试,本文所使用的程序如下:(hello.cs)

 

namespace Hello

{

       class Hello

       {

              public static void Main(string[] args)

              {

                     System.Console.WriteLine("Your name please?");

                     string s = System.Console.ReadLine();

                     Welcome(s);

                     Welcome(s);

              }

 

              public static void Welcome(string name)

              {

                     System.Console.WriteLine("Hello " + name);

              }

 

       }

}

打开命令提示符,进入sscli20目录,键入:

env dbg

进入Rotor的调试环境,如果你还没有BuildRotor的一个Debug版本,那么请参照本系列的第一篇文章来设置你的环境并Build出一个调试版本的Rotor。文章的链接在这里:

.NET Rotor源码研究1 – Building Rotorhttp://blog.youkuaiyun.com/ATField/archive/2006/12/31/1471465.aspx

如果已经Build出来了一个Rotorx86调试版本,那么可以开始动手编译hello.cs (假定hello.cs位于binaries.x86dbg.rotor目录下):

cd binaries.x86dbg.rotor

csc hello.cs

编译之后,启动调试器。这里我们不能直接调试hello.exe,否则.NET将会执行hello.exe,这里我们需要使用clix.exe来运行hello.exe,这样才可以让Rotor来运行hello.exe:

windbg clix hello.exe

请保证WinDbg已经被安装并且在其路径在Path变量中。

程序的加载

启动调试器,我们停在程序加载的位置,Call Stack如下(如果你没有Windows系统DLL所对应的Symbol,那么你看到的会有所不同,这里因为有Symbol,结果更加准确):

ntdll!DbgBreakPoint

ntdll!LdrpDoDebuggerBreak+0x31

ntdll!LdrpInitializeProcess+0xffc

ntdll!_LdrpInitialize+0xf5

ntdll!LdrInitializeThunk+0x10

在本系列的第二篇文章中曾经提到,用到PAL的程序的main实际是在PAL_startup_main,如果你还没有看到第二篇文章的话,连接在这里:

.NET Rotor源码研究2 - PAL http://blog.youkuaiyun.com/ATField/archive/2007/01/12/1481538.aspx

在调试器中输入:

 

bp clix!PAL_startup_main

g

第一条语句的作用是设置断点于clix.exePAL_startup_main函数,第二条语句命令WinDbg继续执行。执行g之后WinDbg很快在clixmain函数停下来,这里的main实际上就是PAL_startup_main,被#define过:

int __cdecl main(int argc, char **argv)

{

// 省略

nExitCode = Launch(pModuleName, pActualCmdLine);

}

 

DWORD Launch(WCHAR* pFileName, WCHAR* pCmdLine)

{

    // 省略

    nExitCode = _CorExeMain2(NULL, 0, pFileName, NULL, pCmdLine);

 

    return nExitCode;

}

这里有不少无关的代码,大部分是分析命令行,直接来到Launch函数调用,Launch函数负责启动ModuleName,也就是hello.exe,启动工作由_CorExeMain2执行。在WindbgF10F11仍然可以工作(当然命令行也可以)。一路执行到_CorExeMain2然后F11,会发现来到了sscoree_CorExeMain2函数,位于sscoree_shims.h之中:

 

SSCOREE_SHIM_RET (

                  __int32,

                  STDMANGLE(_CorExeMain2,20),

                  ( PBYTE   pUnmappedPE,

                    DWORD   cUnmappedPE,

                    LPWSTR  pImageNameIn,

                    LPWSTR  pLoadersFileName,

                    LPWSTR  pCmdLine),

                  ( pUnmappedPE,

                    cUnmappedPE,

                    pImageNameIn,

                    pLoadersFileName,

                    pCmdLine),

                  -1)

这个函数代码很奇怪,只是一些函数调用。仔细观察一下这个头文件,发现这个文件是很有规律的由下面内容组成:

 

SSCOREE_LIB_START (mscorwks)

 

SSCOREE_SHIM_RET (

                  HRESULT,

                  STDMANGLE(MetaDataGetDispenser,12),SSCOREE_LIB_END (mscorwks)

SSCOREE_LIB_END (mscorwks)

 

SSCOREE_LIB_START (mscorpe)

SSCOREE_LIB_END (mscorpe)

 

 

SSCOREE_LIB_START (mscordbi)

这个提示我们SSCOREE.dll会负责将列表中的函数转发到对应的DLL中的对应函数。实际上,这正是sscoree.dll所起到的作用之一,确定Rotor版本,加载对应版本的Rotor,并调用对应版本的Rotor的相应函数,因此sscoree(在.NET中则是mscoree)又被称为Shim。这个SSCOREE_SHIM_RET只是一个宏定义,如下:

#define SSCOREE_SHIM_BODY(FUNC,RET_COMMAND,SIG_RET,SIG_ARGS,ARGS)       /

do {                                                                    /

    SSCOREE_SHIM_CUSTOM_INIT                                            /

    FARPROC proc_addr = SscoreeShimGetProcAddress (                     /

                        SHIMSYM_ ## FUNC,                               /

                        #FUNC);                                         /

    _ASSERTE (proc_addr);                                               /

    if (proc_addr) {                                                    /

        RET_COMMAND ((SIG_RET (STDMETHODCALLTYPE *)SIG_ARGS)proc_addr)ARGS; /

    }                                                                   /

} while (0)

 

 

#define SSCOREE_SHIM_RET(SIG_RET,FUNC,SIG_ARGS,ARGS,ONERROR)            /

extern "C"                                                              /

SIG_RET STDMETHODCALLTYPE FUNC SIG_ARGS                                 /

{                                                                       /

    SSCOREE_SHIM_BODY (FUNC, return, SIG_RET, SIG_ARGS, ARGS);          /

    return ONERROR;                                                     /

}

 

#define SSCOREE_SHIM_NORET(FUNC,SIG_ARGS,ARGS)                          /

extern "C"                                                              /

void STDMETHODCALLTYPE FUNC SIG_ARGS                                    /

{                                                                       /

    SSCOREE_SHIM_BODY (FUNC, ; ,void, SIG_ARGS, ARGS);                  /

}     

可以看到在sscoree中每个类似_CorExeMain2的函数大致作的事情都很类似,首先调用SscoreeShimGetProcAddress获得在Rotor核心DLL中的地址,然后调用之。

回到调试器,按下F11,直接进入SscoreeShimGetProcAddress函数:

FARPROC

SscoreeShimGetProcAddress (

    ShimmedSym SymIndex,

    LPCSTR     SymName)

{

    FARPROC proc;

 

#ifdef TRACE_LOADS

    printf ("SscoreeShimGetProcAddress: Loading Symbol %d (%s)/n",

            SymIndex, g_Syms[SymIndex].Name);

#endif

 

    _ASSERTE (SYM_INDEX_VALID (SymIndex));

    _ASSERTE (SymName);

    _ASSERTE (g_Syms[SymIndex].Name);

    _ASSERTE (!strcmp (g_Syms[SymIndex].Name, SymName));

 

    proc = g_Syms[SymIndex].Proc;

 

    if (proc == NULL) {

        proc = SetupProc(SymIndex, SymName);

    }

 

    return proc;

}

g_Syms是一个全局的数组,用于保存每个函数的实际地址,如果地址=NULL,说明还没有获得此函数的地址,需要调用SetupProc

static

FARPROC

SetupProc (

    ShimmedSym SymIndex,

    LPCSTR     SymName)

{

    HMODULE lib_handle;

    FARPROC proc;

 

    ShimmedLib LibIndex = FindSymbolsLib (SymIndex);

    _ASSERTE (LIB_INDEX_VALID (LibIndex));

 

#ifdef TRACE_LOADS

    printf ("SscoreeShimGetProcAddress: Loading library %d (%S)/n",

            LibIndex, g_Libs[LibIndex].Name);

#endif

 

    lib_handle = g_Libs[LibIndex].Handle;

    if (lib_handle == NULL) {

        lib_handle = SetupLib (LibIndex);

        if (lib_handle == NULL)

            return NULL;

    }

    _ASSERTE (lib_handle);

 

    proc = g_Syms[SymIndex].Proc;

    if (proc == NULL) {

        proc = GetProcAddress (lib_handle, SymName);

        if (!proc) {

#ifdef _DEBUG

            fprintf (stderr,

                        "SscoreeShimGetProcAddress: GetProcAddress (/"%s/") failed/n",

                        SymName);

#endif

            return proc;

        }

        g_Syms[SymIndex].Proc = proc;

    }

 

    return proc;

}

FindSymbols负责找到函数和DLL之间的对应关系:

ShimmedLib

FindSymbolsLib (

    ShimmedSym SymIndex)

{

    // some trickery to figure out which library this symbol is in

    _ASSERTE (SYM_INDEX_VALID (SymIndex));

       

#define SSCOREE_LIB_START(LIBNAME)                                      /

    if (SymIndex < SHIMLIB_ ## LIBNAME) {                               /

        return LIB_ ## LIBNAME;                                         /

    }                                                                   /

    if (SymIndex == SHIMLIB_ ## LIBNAME) {                              /

        return LIB_ ## MAX_LIB;                                         /

    }

#include "sscoree_shims.h"

   

    return LIB_MAX_LIB;

}

这个函数的实现非常有意思,直接定义了两个宏然后includesscoree_shims.h。实际上这是一个很有意思的技巧,sscoree_shims.h中以宏的形式保存了每个函数和每个DLL,这样,通过定义宏的内容,可以对同样的sscoree_shims.h中的内容转换成不同的代码,比如这里就是把这个文件转换成了一系列的if语句,判断函数Index的范围,返回DLL(这里称之为LIB)的Index,避免了重复代码。

再回到SetupProc函数,这次需要注意的SetupProc在调用FindSymbolsLib之后接着调用了SetupLib函数:

 

HMODULE

SetupLib (

    ShimmedLib LibIndex)

{

    HMODULE lib_handle;

    WCHAR FullPath[_MAX_PATH];

 

    if (!PAL_GetPALDirectory (FullPath, _MAX_PATH)) {

        return NULL;

    }

    if (wcslen(FullPath) + wcslen(g_Libs[LibIndex].Name) >= _MAX_PATH) {

        SetLastError(ERROR_FILENAME_EXCED_RANGE);

        return NULL;

    }

    wcsncat(FullPath, g_Libs[LibIndex].Name, _MAX_PATH);

 

    lib_handle = LoadLibrary (FullPath);

    if (lib_handle == NULL) {

#ifdef _DEBUG

        fprintf (stderr,

                    "SscoreeShimGetProcAddress: LoadLibrary (/"%S/") failed/n",

                    FullPath);

        DisplayMessageFromSystem(GetLastError());

#endif

        return lib_handle;

    }

    g_Libs[LibIndex].Handle = lib_handle;

 

#ifdef _DEBUG

    // first time we've hit this library. Run some tests.

    SscoreeVerifyLibrary (LibIndex);

#endif

 

    ROTOR_PAL_CTOR_TEST_RUN(SSCOREE_INT);

 

    return lib_handle;

}

这个函数不长,根据PAL所在目录加载对应的DLL从而实现不同版本的Rotor共存的功能,并且返回加载的DLLHandle。这里我们所需要的DLLmscorwks.dll,是.NET / Rotor 虚拟机的工作站(WorkStation)版本的核心DLL

执行到wcsncat语句之后,在调试器中输入:

 

dv FullPath

这条命令作用是显示FullPath局部变量的值,结果为:

      FullPath = wchar_t [260] "D:/usr/src/sscli20/binaries.x86dbg.rotor/mscorwks.dll"

可以看到我们需要运行mscorwks!_CorExeMain

再度回到SetupProc,这次SetupProc调用GetProcAddress获得对应函数的地址并保存,然后返回。下面的代码就不需要继续执行了。在调试器中输入下面语句:

 

g mscorwks!_CorExeMain2

这条语句让WinDbg执行程序直到遇见mscorwks!_CorExeMain2函数为止:

//*****************************************************************************

// This entry point is called from the native entry piont of the loaded

// executable image.  The command line arguments and other entry point data

// will be gathered here.  The entry point for the user image will be found

// and handled accordingly.

//*****************************************************************************

__int32 STDMETHODCALLTYPE _CorExeMain2( // Executable exit code.

    PBYTE   pUnmappedPE,                // -> memory mapped code

    DWORD   cUnmappedPE,                // Size of memory mapped code

    __in LPWSTR  pImageNameIn,          // -> Executable Name

    __in LPWSTR  pLoadersFileName,      // -> Loaders Name

    __in LPWSTR  pCmdLine)              // -> Command Line

{

 

加载SOS,设置断点

对了,现在我们还需要加载SOS,因为SOS需要mscorwks,因此在这个时候加载SOS正合适。在调试器中输入:

 

.loadby sos mscorwks

这条语句负责将和mscorwks在同一目录下的sos.dll作为WinDbgExtension加载。如果你没有看到任何提示信息,那么加载成功了。如果提示出错,请检查在binaries.x86dbg.rotor目录下面确实存在SOS.dll,并且WinDbg已经被修改过或者MSVCR80D.dll在路径中,具体可以参考本系列第3篇文章:

.NET Rotor源码研究3 - 调试Rotor托管代码的利器:WinDbgSOShttp://blog.youkuaiyun.com/ATField/archive/2007/05/12/1606151.aspx

成功加载之后,为了验证之前我们对IsDebuggerPresent的修改确实生效,输入:

!bpmd hello.exe Hello.Hello.Main

g

前一条命令是SOS命令,负责对Hello.Hello.Main函数设置断点,实际上在CLR中设置断点要比一般程序中设置断点要复杂的多,并且需要notification才可以工作,在后面我将会讲到具体的过程。后面的g命令告诉WinDbg继续执行,注意在WinDbg的输出有如下内容:

(11fc.1088): CLR notification exception - code e0444143 (first chance)

CLR notification: module 'sorttbls.nlp' loaded

(11fc.1088): CLR notification exception - code e0444143 (first chance)

(11fc.1088): CLR notification exception - code e0444143 (first chance)

CLR notification: method 'Hello.Hello.Main(System.String[])' code generated

(11fc.1088): CLR notification exception - code e0444143 (first chance)

JITTED hello!Hello.Hello.Main(System.String[])

若干CLR Notification已经发出,最重要的是最后一个notification,通知Hello.Hello.Main已经被JIT编译成功,之后很快WinDbgHello.Hello.Main函数停下了,说明断点设置成功。

OK,至此我们在Windbg中完整地跟踪了CLIX的启动过程和SSCOREE的函数转发,并成功加载了SOS。下一篇文章将回归mscorwks!_CorExeMain2函数,继续我们在.NET / Rotor之中的探索过程,请继续关注。

 

作者:      张羿/ATField
Blog:     
http://blog.youkuaiyun.com/atfield
转载请注明出处

 

 

% 二叶桨动特性分析(双桨叶五剖面有限元模型),考虑离心力影响 % 旋翼总直径:0.81534m,桨毂直径:50.2mm,自转转速:3800RPM clear; clc; close all; %% 1. 基本参数定义 % 几何参数 rotor_diameter = 0.81534; % 旋翼总直径(m,从一侧桨尖到另一侧) hub_diameter = 0.0502; % 桨毂直径(m) hub_radius = hub_diameter / 2; % 桨毂半径(m) % 单桨叶实际长度 = (旋翼直径 - 桨毂直径) / 2 L_blade = (rotor_diameter - hub_diameter) / 2; n_elem_per_blade = 5; % 每片桨叶单元数量 n_node_per_blade = n_elem_per_blade + 1; % 每片桨叶节点数量 n_blades = 2; % 二叶桨 elem_len = L_blade / n_elem_per_blade; % 每个单元长度(m) % 旋转参数 RPM = 3800; % 自转转速 (RPM) Omega = RPM * 2 * pi / 60; % 角速度 (rad/s) fprintf('旋转角速度: %.2f rad/s (%.0f RPM)\n', Omega, RPM); % 节点编号方案(总节点数11): % - 节点1:公共桨毂中心 % - 桨叶1:节点1, 2, 3, 4, 5, 6 % - 桨叶2:节点1, 7, 8, 9, 10, 11 total_nodes = 1 + (n_node_per_blade - 1) * n_blades; % 存储两片桨叶的节点位置坐标(从桨毂边缘开始) node_pos = cell(n_blades, 1); for b = 1:n_blades % 节点位置从桨毂边缘(hub_radius)到桨尖(hub_radius + L_blade) node_pos{b} = linspace(hub_radius, hub_radius + L_blade, n_node_per_blade); end % 显示几何参数 disp('=== 几何参数 ==='); fprintf('旋翼总直径: %.5f m\n', rotor_diameter); fprintf('桨毂直径: %.5f m\n', hub_diameter); fprintf('单桨叶实际长度: %.5f m\n', L_blade); fprintf('桨叶起始位置(距中心): %.5f m\n', hub_radius); fprintf('桨尖位置(距中心): %.5f m\n', hub_radius + L_blade); fprintf('单元长度: %.5f m\n', elem_len); %% 2. 导入各剖面质量矩阵(6x6,DYMORE convention) M_profiles = cell(n_elem_per_blade, 1); % 单元1(80mm处剖面 - 从桨毂边缘算起) M_profiles{1} = [ 0.289079e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.173399e-13, -0.853452e-03; 0.000000e+00, 0.289079e+00, 0.000000e+00, -0.173399e-13, 0.000000e+00, 0.000000e+00; 0.000000e+00, 0.000000e+00, 0.289079e+00, 0.853452e-03, 0.000000e+00, 0.000000e+00; 0.000000e+00, -0.173399e-13, 0.853452e-03, 0.533859e-04, 0.000000e+00, 0.000000e+00; 0.173399e-13, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.461949e-06, -0.585358e-15; -0.853452e-03, 0.000000e+00, 0.000000e+00, 0.000000e+00, -0.585358e-15, 0.529240e-04; ]; % 单元2(160mm处剖面) M_profiles{2} = [ 0.530206e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, -0.390687e-18, -0.211993e-02; 0.000000e+00, 0.530206e+00, 0.000000e+00, 0.390687e-18, 0.000000e+00, 0.000000e+00; 0.000000e+00, 0.000000e+00, 0.530206e+00, 0.211993e-02, 0.000000e+00, 0.000000e+00; 0.000000e+00, 0.390687e-18, 0.211993e-02, 0.179590e-03, 0.000000e+00, 0.000000e+00; -0.390687e-18, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.155400e-05, 0.222282e-19; -0.211993e-02, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.222282e-19, 0.178036e-03; ]; % 单元3(240mm处剖面) M_profiles{3} = [ 0.407341e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, -0.480828e-13, -0.142755e-02; 0.000000e+00, 0.407341e+00, 0.000000e+00, 0.480828e-13, 0.000000e+00, 0.000000e+00; 0.000000e+00, 0.000000e+00, 0.407341e+00, 0.142755e-02, 0.000000e+00, 0.000000e+00; 0.000000e+00, 0.480828e-13, 0.142755e-02, 0.106001e-03, 0.000000e+00, 0.000000e+00; -0.480828e-13, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.917227e-06, 0.197552e-14; -0.142755e-02, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.197552e-14, 0.105084e-03; ]; % 单元4(320mm处剖面) M_profiles{4} = [ 0.216953e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, -0.655954e-19, -0.554883e-03; 0.000000e+00, 0.216953e+00, 0.000000e+00, 0.655954e-19, 0.000000e+00, 0.000000e+00; 0.000000e+00, 0.000000e+00, 0.216953e+00, 0.554883e-03, 0.000000e+00, 0.000000e+00; 0.000000e+00, 0.655954e-19, 0.554883e-03, 0.300693e-04, 0.000000e+00, 0.000000e+00; -0.655954e-19, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.260190e-06, 0.528450e-20; -0.554883e-03, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.528450e-20, 0.298091e-04; ]; % 单元5(400mm处剖面) M_profiles{5} = [ 0.465444e-01, 0.000000e+00, 0.000000e+00, 0.000000e+00, -0.898566e-20, -0.551386e-04; 0.000000e+00, 0.465444e-01, 0.000000e+00, 0.898566e-20, 0.000000e+00, 0.000000e+00; 0.000000e+00, 0.000000e+00, 0.465444e-01, 0.551386e-04, 0.000000e+00, 0.000000e+00; 0.000000e+00, 0.898566e-20, 0.551386e-04, 0.138398e-05, 0.000000e+00, 0.000000e+00; -0.898566e-20, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.119756e-07, 0.204937e-21; -0.551386e-04, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.204937e-21, 0.137200e-05; ]; %% 3. 导入各剖面刚度矩阵(6x6,Timoshenko梁模型) K_profiles = cell(n_elem_per_blade, 1); % 单元1(80mm处剖面) K_profiles{1} = [ 0.38543891E+08, -0.66963691E-20, -0.56273690E-23, 0.61213913E-25, -0.31799568E-05, -0.11379358E+06; -0.66963691E-20, 0.12953125E+08, -0.21620513E+00, 0.94822720E-03, 0.60167247E-25, 0.54036883E-22; -0.56273690E-23, -0.21620513E+00, 0.24713367E+07, -0.44623520E+04, 0.28654185E-24, -0.14632065E-24; 0.61213913E-25, 0.94822720E-03, -0.44623520E+04, 0.10072960E+03, 0.42614533E-25, -0.24468768E-26; -0.31799568E-05, 0.60167247E-25, 0.28654185E-24, 0.42614533E-25, 0.61593281E+02, 0.33417005E-05; -0.11379358E+06, 0.54036883E-22, -0.14632065E-24, -0.24468768E-26, 0.33417005E-05, 0.70565324E+04; ]; % 单元2(160mm处剖面) K_profiles{2} = [ 0.70694289E+08, -0.87920645E-20, -0.33253495E-23, 0.50586537E-25, -0.55380916E-07, -0.28265461E+06; -0.87920645E-20, 0.23757660E+08, -0.30709474E+00, 0.45802332E-02, -0.39530577E-25, 0.74789743E-22; -0.33253495E-23, -0.30709474E+00, 0.45327500E+07, -0.11084295E+05, 0.72548387E-24, 0.65879888E-26; 0.50586537E-25, 0.45802332E-02, -0.11084295E+05, 0.33885503E+03, 0.11432801E-24, 0.15012058E-27; -0.55380916E-07, -0.39530577E-25, 0.72548387E-24, 0.11432801E-24, 0.20720002E+03, 0.73770970E-05; -0.28265461E+06, 0.74789743E-22, 0.65879888E-26, 0.15012058E-27, 0.73770970E-05, 0.23738259E+05; ]; % 单元3(240mm处剖面) K_profiles{3} = [ 0.54312139E+08, -0.83696149E-20, 0.62221045E-23, 0.63188965E-25, 0.18737318E-04, -0.19033989E+06; -0.83696149E-20, 0.18252228E+08, 0.48802966E-02, -0.68435025E-03, 0.14045370E-25, 0.84467105E-22; 0.62221045E-23, 0.48802966E-02, 0.34823399E+07, -0.74640662E+04, 0.48828173E-24, -0.46403406E-24; 0.63188965E-25, -0.68435025E-03, -0.74640662E+04, 0.20000448E+03, 0.74856810E-25, -0.49685200E-26; 0.18737318E-04, 0.14045370E-25, 0.48828173E-24, 0.74856810E-25, 0.12229702E+03, 0.32527444E-05; -0.19033989E+06, 0.84467105E-22, -0.46403406E-24, -0.49685200E-26, 0.32527444E-05, 0.14011158E+05; ]; % 单元4(320mm处剖面) K_profiles{4} = [ 0.28927020E+08, -0.55019423E-20, -0.19951886E-23, -0.31680102E-25, 0.20243815E-07, -0.73984363E+05; -0.55019423E-20, 0.97212628E+07, 0.97022513E-01, -0.92168379E-04, 0.13735065E-25, 0.37530444E-22; -0.19951886E-23, 0.97022513E-01, 0.18547260E+07, -0.29012603E+04, 0.17499720E-24, 0.24725797E-25; -0.31680102E-25, -0.92168379E-04, -0.29012603E+04, 0.56735317E+02, 0.26440805E-25, 0.51927984E-27; 0.20243815E-07, 0.13735065E-25, 0.17499720E-24, 0.26440805E-25, 0.34692013E+02, -0.31423757E-06; -0.73984363E+05, 0.37530444E-22, 0.24725797E-25, 0.51927984E-27, -0.31423757E-06, 0.39745468E+04; ]; % 单元5(400mm处剖面) K_profiles{5} = [ 0.62059240E+07, -0.11241305E-20, 0.12818189E-24, -0.17512347E-27, 0.53922014E-07, -0.73518077E+04; -0.11241305E-20, 0.20855732E+07, -0.13819523E-02, 0.24279922E-06, 0.42246092E-27, 0.33534716E-23; 0.12818189E-24, -0.13819523E-02, 0.39790589E+06, -0.28829546E+03, 0.13609137E-25, 0.71224377E-27; -0.17512347E-27, 0.24279922E-06, -0.28829546E+03, 0.26113084E+01, 0.12534636E-26, 0.17057581E-29; 0.53922014E-07, 0.42246092E-27, 0.13609137E-25, 0.12534636E-26, 0.15967422E+01, 0.40783515E-08; -0.73518077E+04, 0.33534716E-23, 0.71224377E-27, 0.17057581E-29, 0.40783515E-08, 0.18293344E+03; ]; %% 4. 计算离心力引起的几何刚度矩阵(初应力刚度矩阵) K_geo_profiles = cell(n_elem_per_blade, 1); for i = 1:n_elem_per_blade % 从质量矩阵提取单元质量信息 mass_elem = M_profiles{i}(1,1); % 轴向质量 Iy_elem = M_profiles{i}(5,5); % 绕y轴转动惯量 Iz_elem = M_profiles{i}(6,6); % 绕z轴转动惯量 % 单元位置(距旋转中心)- 取单元中点位置 nodeA_pos = node_pos{1}(i); nodeB_pos = node_pos{1}(i+1); r_mid = (nodeA_pos + nodeB_pos) / 2; % 单元中点半径 % 计算离心力: F = m * r * Ω² F_centrifugal = mass_elem * r_mid * Omega^2; % 构建Timoshenko梁单元的几何刚度矩阵(考虑离心力) % 参考旋转梁理论推导的几何刚度矩阵形式 K_geo_elem = zeros(6,6); % 轴向刚度增强(由离心力引起) K_geo_elem(1,1) = F_centrifugal / elem_len; % 弯曲刚度增强(考虑转动惯量影响) K_geo_elem(5,5) = (Iy_elem * Omega^2) / elem_len; % 绕y轴 K_geo_elem(6,6) = (Iz_elem * Omega^2) / elem_len; % 绕z轴 % 剪切刚度影响 K_geo_elem(2,2) = 0.1 * F_centrifugal / elem_len; % 经验系数 K_geo_elem(3,3) = 0.1 * F_centrifugal / elem_len; % 经验系数 % 扭转刚度影响 K_geo_elem(4,4) = 0.05 * F_centrifugal * r_mid / elem_len; % 经验系数 K_geo_profiles{i} = K_geo_elem; end % 显示刚度矩阵特性 disp('=== 刚度矩阵特性 ==='); for i = 1:n_elem_per_blade fprintf('单元%d - 结构刚度最大项: %.3e, 几何刚度最大项: %.3e\n', ... i, max(max(abs(K_profiles{i}))), max(max(abs(K_geo_profiles{i})))); end %% 5. 组装整体质量矩阵、刚度矩阵和几何刚度矩阵 DOF_per_node = 6; % 每个节点6个自由度 DOF_total = DOF_per_node * total_nodes; % 总自由度:11节点×6=66 M_global = zeros(DOF_total); % 整体质量矩阵 K_global = zeros(DOF_total); % 整体刚度矩阵 K_geo_global = zeros(DOF_total); % 整体几何刚度矩阵 % 为每片桨叶定义完整节点列表(包含根部) blade_nodes = { [1, 2, 3, 4, 5, 6]; % 桨叶1节点 [1, 7, 8, 9, 10, 11]; % 桨叶2节点 }; % 循环组装两片桨叶的单元 for b = 1:n_blades % 桨叶编号:1(左)、2(右) for i = 1:n_elem_per_blade % 每个桨叶的5个单元 % 从节点列表中获取当前单元的两个节点 nodeA = blade_nodes{b}(i); nodeB = blade_nodes{b}(i+1); % 计算自由度索引 dofA = (nodeA - 1)*DOF_per_node + 1 : nodeA*DOF_per_node; dofB = (nodeB - 1)*DOF_per_node + 1 : nodeB*DOF_per_node; % 检查索引范围 if any(dofA > DOF_total) || any(dofB > DOF_total) error('自由度索引超出范围: 桨叶=%d, 单元=%d, nodeA=%d, nodeB=%d', b, i, nodeA, nodeB); end % 获取单元矩阵 M_elem = M_profiles{i}; K_elem = K_profiles{i} / elem_len; % 结构刚度矩阵 K_geo_elem = K_geo_profiles{i}; % 几何刚度矩阵 % 组装质量矩阵 M_global(dofA, dofA) = M_global(dofA, dofA) + M_elem; M_global(dofB, dofB) = M_global(dofB, dofB) + M_elem; % 组装结构刚度矩阵 K_global(dofA, dofA) = K_global(dofA, dofA) + K_elem; K_global(dofB, dofB) = K_global(dofB, dofB) + K_elem; K_global(dofA, dofB) = K_global(dofA, dofB) - K_elem; K_global(dofB, dofA) = K_global(dofB, dofA) - K_elem; % 组装几何刚度矩阵 K_geo_global(dofA, dofA) = K_geo_global(dofA, dofA) + K_geo_elem; K_geo_global(dofB, dofB) = K_geo_global(dofB, dofB) + K_geo_elem; K_geo_global(dofA, dofB) = K_geo_global(dofA, dofB) - K_geo_elem; K_geo_global(dofB, dofA) = K_geo_global(dofB, dofA) - K_geo_elem; end end %% 6. 应用边界条件(桨毂中心固定) fixed_dof = 1:DOF_per_node; % 节点1(桨毂中心)的6个自由度全部固定 free_dof = setdiff(1:DOF_total, fixed_dof); % 自由自由度 % 缩减矩阵 K_red = K_global(free_dof, free_dof); M_red = M_global(free_dof, free_dof); K_geo_red = K_geo_global(free_dof, free_dof); % 组合刚度矩阵(结构刚度 + 几何刚度) K_combined_red = K_red + K_geo_red; %% 7. 求解固有频率和振型(考虑与不考虑离心力两种情况) % 情况1: 不考虑离心力 [V_no_cent, D_no_cent] = eig(K_red, M_red); lambda_no_cent = diag(D_no_cent); frequencies_no_cent = sqrt(abs(lambda_no_cent)) / (2*pi); [frequencies_no_cent, idx_no_cent] = sort(frequencies_no_cent); V_no_cent = V_no_cent(:, idx_no_cent); % 情况2: 考虑离心力 [V_with_cent, D_with_cent] = eig(K_combined_red, M_red); lambda_with_cent = diag(D_with_cent); frequencies_with_cent = sqrt(abs(lambda_with_cent)) / (2*pi); [frequencies_with_cent, idx_with_cent] = sort(frequencies_with_cent); V_with_cent = V_with_cent(:, idx_with_cent); %% 8. 结果输出与对比 disp(' '); disp('=== 模态分析结果对比 ==='); disp('(考虑与不考虑3800RPM离心力作用)'); fprintf('旋翼总直径: %.5f m\n', rotor_diameter); fprintf('单桨叶实际长度: %.5f m\n', L_blade); disp(' '); % 创建对比表格 disp('阶数 | 不考虑离心力 (Hz) | 考虑离心力 (Hz) | 变化量 (Hz) | 变化率 (%)'); disp('---------------------------------------------------------------------'); for i = 1:min(10, length(frequencies_no_cent)) freq1 = frequencies_no_cent(i); freq2 = frequencies_with_cent(i); delta = freq2 - freq1; delta_percent = (delta / freq1) * 100; fprintf('%4d | %19.4f | %18.4f | %11.4f | %10.2f\n', ... i, freq1, freq2, delta, delta_percent); end %% 9. 振型可视化(对比两种情况) dof_indices = [2, 3, 4, 5, 6]; % 2=剪切1, 3=剪切2, 4=扭转, 5=弯曲1, 6=弯曲2 dof_name = {'剪切方向1', '剪切方向2', '扭转方向', '弯曲方向1', '弯曲方向2'}; for fig = 1:length(dof_indices) % 创建对比图 figure('Name', ['振型对比 - ', dof_name{fig}], 'Position', [100, 100, 1000, 800]); for mode = 1:3 % 前3阶模态 % 不考虑离心力的振型 subplot(3, 2, 2*mode-1); % 提取左桨叶(桨叶1)振型 disp_left = zeros(1, n_node_per_blade); pos_left = zeros(1, n_node_per_blade); for node_idx = 1:n_node_per_blade actual_node = blade_nodes{1}(node_idx); dof_idx = (actual_node - 1)*DOF_per_node + dof_indices(fig); pos_left(node_idx) = node_pos{1}(node_idx); if ismember(dof_idx, free_dof) red_idx = find(free_dof == dof_idx); disp_left(node_idx) = V_no_cent(red_idx, mode); else disp_left(node_idx) = 0; end end % 提取右桨叶(桨叶2)振型 disp_right = zeros(1, n_node_per_blade); pos_right = zeros(1, n_node_per_blade); for node_idx = 1:n_node_per_blade actual_node = blade_nodes{2}(node_idx); dof_idx = (actual_node - 1)*DOF_per_node + dof_indices(fig); pos_right(node_idx) = node_pos{2}(node_idx); if ismember(dof_idx, free_dof) red_idx = find(free_dof == dof_idx); disp_right(node_idx) = V_no_cent(red_idx, mode); else disp_right(node_idx) = 0; end end % 归一化 max_amp = max([abs(disp_left), abs(disp_right)]); if max_amp > 0 disp_left = disp_left / max_amp; disp_right = disp_right / max_amp; end % 绘图 rectangle('Position', [-hub_radius, -1.2, hub_diameter, 2.4], ... 'FaceColor', [0.8 0.8 0.8], 'EdgeColor', 'k', 'LineWidth', 1); hold on; plot(pos_left, disp_left, 'o-', 'LineWidth', 1.5, 'Color', [0.2 0.5 0.8], 'DisplayName', '左桨叶'); plot(-pos_right, disp_right, 's-', 'LineWidth', 1.5, 'Color', [0.8 0.3 0.3], 'DisplayName', '右桨叶'); plot(0, 0, 'ko', 'MarkerSize', 8, 'MarkerFaceColor', 'k'); xlabel('距中心距离(m)'); ylabel('归一化振幅'); title(['第', num2str(mode), '阶 (无离心力), f=', num2str(frequencies_no_cent(mode), '%.4f'), 'Hz']); legend('Location', 'best'); grid on; x_range = max(hub_radius + L_blade) * 1.1; xlim([-x_range, x_range]); ylim([-1.2, 1.2]); plot([hub_radius, hub_radius], [-1.2, 1.2], 'k--', 'LineWidth', 0.5); plot([-hub_radius, -hub_radius], [-1.2, 1.2], 'k--', 'LineWidth', 0.5); box on; hold off; % 考虑离心力的振型 subplot(3, 2, 2*mode); % 提取左桨叶(桨叶1)振型 disp_left_cent = zeros(1, n_node_per_blade); for node_idx = 1:n_node_per_blade actual_node = blade_nodes{1}(node_idx); dof_idx = (actual_node - 1)*DOF_per_node + dof_indices(fig); if ismember(dof_idx, free_dof) red_idx = find(free_dof == dof_idx); disp_left_cent(node_idx) = V_with_cent(red_idx, mode); else disp_left_cent(node_idx) = 0; end end % 提取右桨叶(桨叶2)振型 disp_right_cent = zeros(1, n_node_per_blade); for node_idx = 1:n_node_per_blade actual_node = blade_nodes{2}(node_idx); dof_idx = (actual_node - 1)*DOF_per_node + dof_indices(fig); if ismember(dof_idx, free_dof) red_idx = find(free_dof == dof_idx); disp_right_cent(node_idx) = V_with_cent(red_idx, mode); else disp_right_cent(node_idx) = 0; end end % 归一化 max_amp_cent = max([abs(disp_left_cent), abs(disp_right_cent)]); if max_amp_cent > 0 disp_left_cent = disp_left_cent / max_amp_cent; disp_right_cent = disp_right_cent / max_amp_cent; end % 绘图 rectangle('Position', [-hub_radius, -1.2, hub_diameter, 2.4], ... 'FaceColor', [0.8 0.8 0.8], 'EdgeColor', 'k', 'LineWidth', 1); hold on; plot(pos_left, disp_left_cent, 'o-', 'LineWidth', 1.5, 'Color', [0.2 0.5 0.8], 'DisplayName', '左桨叶'); plot(-pos_right, disp_right_cent, 's-', 'LineWidth', 1.5, 'Color', [0.8 0.3 0.3], 'DisplayName', '右桨叶'); plot(0, 0, 'ko', 'MarkerSize', 8, 'MarkerFaceColor', 'k'); xlabel('距中心距离(m)'); ylabel('归一化振幅'); title(['第', num2str(mode), '阶 (有离心力), f=', num2str(frequencies_with_cent(mode), '%.4f'), 'Hz']); legend('Location', 'best'); grid on; x_range = max(hub_radius + L_blade) * 1.1; xlim([-x_range, x_range]); ylim([-1.2, 1.2]); plot([hub_radius, hub_radius], [-1.2, 1.2], 'k--', 'LineWidth', 0.5); plot([-hub_radius, -hub_radius], [-1.2, 1.2], 'k--', 'LineWidth', 0.5); box on; hold off; end sgtitle(['双桨叶振型对比 - ', dof_name{fig}, ' (3800 RPM离心力影响)']); end 帮我分析一下这个代码
最新发布
08-29
评论
成就一亿技术人!
拼手气红包6.0元
还能输入1000个字符
 
红包 添加红包
表情包 插入表情
 条评论被折叠 查看
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值